Predatory fish depletion and recovery potential on Caribbean reefs.

The natural, prehuman abundance of most large predators is unknown because of the lack of historical data and a limited understanding of the natural factors that control their populations. Determining the supportable predator biomass at a given location (that is, the predator carrying capacity) would help managers to optimize protection and would provide site-specific recovery goals. We assess the relationship between predatory reef fish biomass and several anthropogenic and environmental variables at 39 reefs across the Caribbean to (i) estimate their roles determining local predator biomass and (ii) determine site-specific recovery potential if fishing was eliminated. We show that predatory reef fish biomass tends to be higher in marine reserves but is strongly negatively related to human activities, especially coastal development. However, human activities and natural factors, including reef complexity and prey abundance, explain more than 50% of the spatial variation in predator biomass. Comparing site-specific predator carrying capacities to field observations, we infer that current predatory reef fish biomass is 60 to 90% lower than the potential supportable biomass in most sites, even within most marine reserves. We also found that the scope for recovery varies among reefs by at least an order of magnitude. This suggests that we could underestimate unfished biomass at sites that provide ideal conditions for predators or greatly overestimate that of seemingly predator-depleted sites that may have never supported large predator populations because of suboptimal environmental conditions.

INTRODUCTION

Overfishing and habitat degradation have caused the decline of countless large marine predator species from most of the world’s oceans (–). Their widespread depletion has indirectly modified (or eliminated) species interactions, redistributed the flow of energy, and altered ecosystem functioning and services (, , –). Fishing has substantially changed the structure of fish communities in coral reefs, particularly reducing the abundance and size of large predatory fishes (–), including widespread declines of reef sharks (–). However, quantifying the impacts of fishing has proven challenging because we generally lack quantitative, spatially replicated baseline data on the preexploitation state of fish assemblages. Historical analyses suggest extensive reduction of large marine predators that typically preceded population assessments (–). The extent of fishing on coral reefs has been estimated using fishing intensity gradients based on catch efforts (, , , , ), fishing mortality differences between inside and outside marine reserves (, , –), and human pressure indices as a proxy for fishing intensity (, , , , –). In general, these studies show strong spatial and temporal gradients in fishing pressure and fish biomass and declines of roughly 50 to 95% compared to presumed unfished biomass values.

Although returning many predator populations to pristine historical baseline abundance is unrealistic, the knowledge of the current natural state of predator assemblages in the absence of fishing (that is, potential carrying capacity) gives managers reasonable science-based recovery goals to evaluate the efficacy of management. However, ecological baselines vary with environmental context (). Therefore, to assess the degree to which human activities have altered marine communities and to estimate local and regional potential carrying capacity, we need a better understanding of the factors that control the structure and composition of unexploited communities. We know little about the natural abundance and distribution of predator assemblages across landscapes and regions (, , , ). We tend to assume that predators used to be ubiquitous—present at all locations (, )—but our knowledge of their natural spatiotemporal distribution and abundance is limited. Given the influence of environmental characteristics on predator populations, the assumption that all sites and regions have the potential, in the absence of fishing, to sustain fish communities similar to “quasi-pristine” baseline sites (, ) may be unfounded.

Home-range and foraging patterns of marine predators (for example, sharks, tuna, and turtles) respond to prey availability (–) and can also be influenced by other predators, competitors, temperature, habitat structural complexity, and numerous other biotic and abiotic factors (, ). For example, sea surface temperature and primary productivity cause differences in reef shark abundances across coral reefs of the central Pacific (). In the western Pacific, ~19 to 53% of the variation in diversity and biomass of large-bodied reef fish is explained by atoll position, temperature, depth, wave energy, distance to deep water, and topography complexity, whereas fishing explains ~26 to 60% (, ). Yet, in exploited ecosystems, bottom-up forcing can be difficult to detect because predators are affected by spatiotemporal variation in top-down control by humans (that is, hunting or fishing) that obscures any response to natural environmental differences among sites (, , , ).

Here, we quantify the spatial relationships between human population attributes (as proxies for fishing intensity and other anthropogenic impacts), prey abundance and other environmental characteristics, and the biomass of predatory reef fishes (sharks, grouper, etc.) at 39 reefs across the central-western Caribbean. We used this information to determine (i) site-specific carrying capacities and recovery potential, (ii) the degree of predatory fish biomass missing from Caribbean reefs, and (iii) the effectiveness of marine reserves in protecting predatory reef fishes. Our results not only indicate high recovery potential of predatory fish biomass on Caribbean reefs but also suggest that natural predator abundance varies greatly among sites because of environmental characteristics, including resource availability and habitat heterogeneity.

RESULTS

Spatial variation of observed predatory fish biomass

Observed reef fish and predatory fish biomass varied greatly among sites even within marine reserves. Among the 39 fore-reef sites studied, the average total fish biomass per site (mean ± SE) ranged from ~55 ± 8 g m−2 at Mexico Rocks in Belize to ~445 ± 148 g m−2 at El Peruano inside the Gardens of the Queen marine reserve in Cuba (Fig. 1). Mean total fish biomass was ~137 ± 9 g m−2 on unprotected reefs and ~237 ± 16 g m−2 across marine reserves (Fig. 1 and fig. S1). The mean biomass of apex predators and piscivore-invertivores combined (hereafter “predators”) ranged from ~3 ± 1 g m−2 at Ebano to ~260 ± 78 g m−2 at El Peruano, both sites in Cuba (Fig. 1). Although there was high spatial variation among reefs within countries and protection levels, total fish and predator biomass was higher within the marine reserves of Abaco, Cuba, and Mexico than in reserves of Belize (Fig. 1 and fig. S1). The proportion of trophic guilds varied across sites (Fig. 1 and fig. S2), but the biomass of all lower trophic levels was positively correlated with predator biomass (rs = ~0.20 to 0.35, P < 0.001; fig. S3). For detailed description of fish biomass spatial difference among sites, see section S1.

Fig. 1

Biomass of reef fishes by sites and trophic level.

Mean biomass of trophic guilds per reef site + 1 SE for total fish biomass. Sites are organized from low to high total fish biomass. Note that the observed biomass distribution is not meant to be representative of the Caribbean because sites were not randomly selected, that is, the higher biomass sites were specifically targeted to quantify reef fish structure at locations with limited fishing. Moreover, we did not sample the Caribbean regions where overfishing is the most severe, for example, Jamaica. Trophic categories were based on dietary information. No-take zones (NTZs) and minimum fished marine protected areas (MPAs) are noted as reserves (gray dots). For site abbreviations and protection level information, see table S2. For species list in each group, see table S3.

Covariates explaining the spatial variation of observed predatory fish biomass

Predatory fish biomass was strongly and negatively correlated with human-related variables (for example, coastal development and human population density) but positively related to habitat reef complexity (Fig. 2, Table 1, and fig. S4). The combination of top models (ΔAICc < 2, where ΣwAICc > 0.95) analyzing all sites (model set A) and sites within marine reserves (model set B) included “coastal development” for all predator groups (Fig. 2A and Table 1). “Human population” was not selected in the top models across all sites (Table 1); however, it was a good predictor for predator biomass within marine reserves (Fig. 2A). Only the biomass of apex predators was negatively related to “cultivated land” (Fig. 2A and Table 1). Among the physical cofactors, only the biomass of apex predators showed a weak positive relationship with “ocean productivity,” whereas the quadratic term of “temperature” was not important (Fig. 2B, Table 1, and fig. S4). In contrast, the biomass of both predatory fish was positively related to “reef complexity” (Fig. 2B, Table 1, and figs. S4 and S5).

Fig. 2

Coefficient estimates of generalized mixed-effects models.

Mean coefficient estimates (±95% confidence interval) of top models (ΔAICc < 2, where ΣwAICc > 0.95) for apex predators, piscivore-invertivores (Pisc-invertivores), and total predators. Longer confidence intervals are truncated to improve visualization. Black circles and lines are estimates of the models that included all study sites (model set A). Blue squares and lines are for models within marine reserves (model set B). Gray horizontal lines divide variables by anthropogenic (A), physical (B), biological (C), and management categories (D). Only estimates that improved model fit are shown. Poaching high, high level of poaching.

Table 1

Summary of generalized linear mixed-effects models.

Results of the generalized linear mixed-effects model (GLMM) comparisons using Akaike’s information criterion corrected for small sample sizes (AICc) for apex predators, piscivore-invertivores, and total predators. Only the null model, an exploratory model that outperformed the null model, and final models (ΔAICc < 2, where ΣwAICc > 0.95) are shown. Model set A included all sites, whereas model set B included sites within reserves. Parameters are as follows: model maximum log-likelihood (LL), degrees of freedom (df), change in AICc (ΔAICc), AICc weights (wAICc), and pseudocoefficient of determination for GLMMs (R2). Models are ordered by increasing wAICc, and the final models show bolded wAICc. Model covariates include the following: Ag, reserve age; Al, macroalgae; Db, distance to reef break; De, depth; Dm, distance to mangrove; Cd, coastal development within 50 km; Cl, cultivated land within 50 km; Co, corals; He, herbivores; Hu, human population density within 50 km; In, invertivores; Go, gorgonians; Ma, mangrove perimeter within 5 km; Tp2, quadratic term of minimum monthly mean sea surface temperature; Om, omnivores; Pi, piscivore-invertivores; Pl, planktivores; Po, poaching level within reserve; Pp, net primary production; Pr, protection level; Ra, reef area within 5 km; Ru, reef complexity; Si, reserve size; We, wave exposure. See table S4 for units.

Models	LL	df	ΔAICc	wAICc	R2
Apex predators
Set A
Null	−324.9	6	20.16	0.00
Cd + Tp2 + Ma + Co + Go + Pi + Pr	−312.4	14	12.00	0.00	0.49
Cd + Ru + Co + Pi + Pr	−309.1	12	1.10	0.18	0.49
Cl + Pp + Ru + Co + Pi + Pr	−307.7	13	0.40	0.25	0.50
Cl + Ru + Co + Pi	−310.9	10	0.35	0.26	0.49
Cd + Tp2 + Ru + Co + Pi	−309.6	11	0.00	0.31	0.49
Set B
Null	−129.5	6	14.3	0.00
Cd + Ru + Ma + Go + Pi + In + Om	−121.3	13	14.0	0.00	0.45
Ru + Ma + Co + Pi + Po	−117.2	11	1.09	0.14	0.42
Hu + Ru + Pi	−119.4	9	0.68	0.17	0.42
Cl + Ru + Co + Pi	−118.0	10	0.28	0.21	0.42
Pp + Ru + Co + Po	−117.9	10	0.12	0.23	0.43
Cd + Ru + Co + Go + Pi	−116.7	11	0.00	0.24	0.44
Piscivore-invertivore
Set A
Null	−252.6	5	65.68	0.00
Cd + Pp + Tp2 + We + Ru + Ma + Co + Al + Go + In + Om + Pl + He + Pr	−207.1	20	6.74	0.03	0.57
Cd + Tp2 + Ru + Co + Go + In + Om + Pl + He + Pr	−208.1	16	0.00	0.97	0.56
Set B
Null	−85.91	5	18.81	0.00
Cd + Tp2 + We + Ru + Co + Go + In + Om + Pl + He + Po	−65.69	17	7.96	0.01	0.58
Hu + Ru + Ma + Co + In + Om + Po	−68.29	12	0.64	0.42	0.54
Cd + Ru + Co + In + Om + Po	−69.85	11	0.00	0.57	0.55
Total predators
Set A
Null	−249.0	5	45.69	0.00
Cd + Pp + Tp2 + We + De + Ru + Db + Dm + Ma + Co + Al + Go + In + Om + Pl + He + Pr	−219.2	21	19.77	0.00	0.60
Cd + Tp2 + Ru + Co + Go + In + Om + Pl + He + Pr	−214.9	16	0.00	1.00	0.61
Set B
Null	−84.29	5	8.14	0.00
Cd + Ru + Ma + Co + Go + In + Om + Pl + He + Si + Ag + Po	−71.49	15	6.70	0.00	0.56
Hu + Ru + Co + In + Om + Po	−73.24	11	2.30	0.23	0.53
Cd + Ru + Co + In + Om + Po	−74.39	11	0.00	0.77	0.55

Predatory fish biomass was more related to potential prey than to habitat features, such as coral and gorgonian cover or connectivity to mangrove (Fig. 2C and Table 1). For example, only the biomass of apex predators tended to be higher at few reefs within marine reserves associated with mangrove (Fig. 2C). The biomass of both predator groups was mostly unrelated to “coral cover” or “gorgonian abundance” (Fig. 2C), although piscivore-invertivores were negatively correlated to the latter across all sites (Fig. 2C and Table 1). In contrast, lower trophic fish groups (except for planktivores) drove some of the spatial variation of predator biomass across sites (Fig. 2C and Table 1). Particularly, the biomass of piscivore-invertivores was positively and strongly correlated with the biomass of invertivore, omnivore, and herbivore fishes (except for herbivores within marine reserves) (Fig. 2C and Table 1).

Both predatory fish groups responded differently to protection level (for example, no protection, MPA, and NTZ) (Fig. 2D and Table 1). No protection was used to set the comparisons for the NTZ and MPA categories. The biomass of apex predators was not related to NTZs (that is, no fishing) and showed no correlation with MPAs (that is, some fishing is allowed) across sites (Fig. 2D). In contrast, the biomass of piscivore-invertivores showed no association with NTZ but a strong negative relation with MPAs (Fig. 2D). Finally, apex predator biomass tended to be lower within marine reserves with higher “poaching levels” (Fig. 2D). Reserve age and size were not important for any predator group in our study (Table 1).

Overall, the models analyzing all sites (model set A) “explained” ~50, ~57, and ~61% of the spatial variation among sites in the biomass of apex predators, piscivore-invertivores, and total predators, respectively (Table 1). Within marine reserves, all these covariates explained ~43% of the variation of apex predators among sites, ~56% of that of piscivore-invertivores, and ~58% of that of total predators (Table 1). For detailed results of the relationship of predatory fish biomass and main cofactors, see section S2.

Reconstructed potential carrying capacity for predatory fish biomass

The potential site-specific carrying capacity range for predatory fish biomass was estimated by comparing the difference between expected (light blue) and observed (orange) mean biomass after setting human-related variables to zero and considering all sites protected (Fig. 3). We inferred that three of four reefs (30 of 39 sites), even within marine reserves, support 60 to 90% less predatory fish biomass than the potential carrying capacity at a given site, likely because of human activities associated with coastal development, such as fishing (Fig. 3; see table S1 for values). For example, we observed that the average biomass of apex predators (for example, sharks, jacks, barracudas, tarpon, and large-bodied groupers) was less than 10% of the estimated potential carrying capacity, whereas the observed average biomass of piscivore-invertivores (for example, medium-bodied fish, such as jacks, snappers, and groupers) was ~18% of the estimated biomass in the absence of fishing across sites.

Fig. 3

Estimated biomass of predatory reef fish in the absence of human activities.

Boxplot of the observed (orange) and predicted (light blue) median (50 and 99% quartiles) of predatory reef fish biomass across survey sites (ordered from lowest to highest biomass). Horizontal dashed lines are the observed (orange) and predicted (light blue) regional means. Predicted biomass was based on the best explanatory model given, with no coastal development within 50 km (that is, in the absence of humans) and every site considered as an NTZ (that is, no fishing). On the basis of the predictive models, 60 to 90% of the potential predatory fish biomass has been extracted in three of four reefs. NTZs and MPAs with minimum fishing are noted as marine reserves (gray dots). For better representation, y axis is in log scale. For site codes, see table S2.

Few sites within marine reserves showed observed biomass averages within the expected range in the absence of human activities (Fig. 3 and table S1). For example, total predatory fish biomass at sites within the marine reserves of Gardens of the Queen in Cuba and Dry Tortugas in Florida showed predatory fish biomass within values predicted by the models (Fig. 3). In contrast, reef sites within other marine reserves, such as Hol Chan and Half Moon in Belize, regarded as well managed, showed substantially lower observed than expected values (Fig. 3 and table S1). On average, sites within NTZs showed a ~66% less than expected predatory fish biomass, whereas there was ~88% less at sites with no protection and within MPAs.

DISCUSSION

Reef complexity and lower trophic fish were strong predictors of predatory fish biomass. Habitat complexity is thought to attract predators, in part, because it is often associated with prey richness and abundance (). We found that the standing biomass of lower trophic levels was positively related to predatory fish biomass, especially for piscivore-invertivores (Fig. 3). Other environmental variables, such as productivity, temperature, and benthic habitat features played a less important role (for a discussion on the rest of variables, see section S3). Our models explained more than 50% of the observed spatial variation in predator biomass (Table 1). However, we caution that there are likely additional important variables that we did not consider. For example, larval supply (), intraguild competition and predation (, ), and habitat connectivity (, ) also regulate predatory fish populations.

Comparing predicted biomass of predatory reef fish in the absence of humans (that is, no fishing) to observed values suggests that fishing and other human impacts may have reduced Caribbean predatory reef fish biomass to 10 to 40% of supportable levels (Fig. 3 and table S1). These differences between observed and predicted biomass may not necessarily represent a decline from historical levels but a departure of observed biomass from potential carrying capacity that is likely due to fishing and other human activities. These marked differences, not previously documented for the entire assemblage of Caribbean predatory reef fishes, although perhaps not unexpected, are concordant with other large-scale studies that show losses of more than 90% from the assumed baselines in oceanic and coastal waters across the globe, primarily due to overfishing (, , ). For example, comparisons of fished reefs with undisturbed reefs across gradients of exploitation suggest total reef fish declines of ~75 to 90% across the Caribbean (, ), ~65 to 93% decline in fished Kenyan reefs in comparison with the underexploited Maldives (, , ), ~70 to 93% decline in the main Hawaiian islands (, ), and ~60 to 90% decline across the central Pacific (, , ). Our model accounts for current differences in environmental conditions among sites that also strongly influenced predatory reef fish assemblages (, , ). Thus, we avoided a common pitfall of using remote locations or well-established fisheries closures with distinct environmental conditions to infer general community carrying capacity. The mean predicted supportable predatory fish biomass ranged by a factor of 70 (26 to 1847 g m−2) among sites, considering both extremes (Fig. 3 and table S1), emphasizing the large degree of context dependency on predatory reef fish carrying capacity.

On the basis of our model predictions, some fore-reef sites may be potential hot spots for predatory fish biomass, with predicted average biomass of more than 500 g m−2, if anthropogenic stressors were eliminated and fishing regulations were better enforced. For example, our model predicts that Columbia Reef within fisheries closures of Cozumel, Mexico, could support, on average, 10 times (~725 g m−2) the current levels of predatory fish biomass (Fig. 3 and table S1). The central and north sites of Banco Chinchorro in Mexico could support average predator fish biomass of more than 800 g m−2. The observed values at these sites were only ~10% of their predicted supportable biomass (Fig. 3 and table S1). Fished sites, such as Bacunayagua in the northern site of Cuba and Rocky Point in the south tip of Abaco, Bahamas, could potentially reach 10 and 5 times higher biomass than the current levels, respectively (Fig. 3 and table S1). This information could be used by managers as benchmarks to trigger conservation initiatives (for example, where to place new marine reserves) and to better tailor conservation efforts for strategic protection and restoration. However, predicted biomass in our study should be interpreted with caution because there are several ecological variables (for example, competition and intraguild predation) and other environmental factors that we did not include in the models that may influence realized site-specific fish biomass.

Predatory fish biomass predicted in our study is also considerably higher than unfished fish biomass from some recent studies, which ranges from 100 to 150 g m−2 across reefs of the Caribbean, Indian, and Pacific oceans (–, ). Several reasons can explain these discrepancies. First, we included large predatory reef fishes in our surveys, which greatly contribute to the total fish biomass. However, most of these studies excluded apex predators, such as sharks and large jacks, and the survey methodology was not designed to account for these predators (–, ). Second, we used noninstantaneous surveys, recording all fishes within, passing into, or across the survey area in front of the divers as they swim along the transect. Evidence from simulation models suggests that noninstantaneous surveys tend to overestimate the abundance of large and vagile fishes (). The methodology of some of these studies includes instantaneous surveys () that are not comparable to our study design. However, our observed values are very similar to studies with comparable goals and methods (). Third, reef fish biomass values estimated in these studies are averages from several sites across a wide range of reef habitats and depths (, , ) and are not necessarily representative of site-specific and structurally complex fore-reef habitats at a depth of 10 to 15 m as in our study. Therefore, we reiterate that our fish abundance estimates (including predicted predatory fish biomass) are site-specific for fore-reef habitats and variable according to their environmental conditions.

Given the sensitivity of predatory reef fish surveys to methodology, observer biases, and environmental conditions, the specific carrying capacity values predicted for our sites should not be applied to other reefs or habitat types. Also, these site-specific supportable biomass estimates are potential values under current conditions. We cannot know precisely the past, prehuman predator biomass of a reef because many of the key environmental drivers have presumably changed. Thus, potential carrying capacities are both time- and location-specific, and they are not necessarily indicators of the past. Our results have conservation applications because they demonstrate how much predatory fish biomass may be missing, how ineffective most marine reserves and MPAs are for both top predators and mesopredators (Fig. 2), and how much more potential for recovery and fisheries productivity still exists (Fig. 3). We cannot predict how much biomass would be achieved by restricting fishing at a given site, but we can say with some certainty which sites are likely to support the largest predatory fish populations. This approach could also be useful in providing site-specific targets in data-poor multispecies reef fisheries, where historical and time series information is absent.

Predatory fishes represented more than 40% of the total fish biomass at some of the relatively isolated reefs and inside the large and well-enforced marine reserves that we surveyed (Fig. 1). Most protected reefs had higher total fish and predator biomass (Fig. 1), and the abundance of apex predators, such as sharks, groupers, snappers, and jacks, declined across a gradient of human effects (Fig. 3). Most of these predators were entirely absent from unprotected sites (Fig. 1), a finding concordant with presence-absence surveys performed by citizen scientists in the region () and scientific surveys (). Large reef fish predators are rare throughout the Caribbean and occupy only a small fraction of sites because of selective targeting by fishermen and greater vulnerability to even low levels of fishing (, ). Moreover, because of their slower growth and lower fecundity rates, large fishes recover more slowly once a reserve is established—often requiring decades to become reestablished (, ). These factors lead to nonrandom changes in composition and reductions in richness as fishing pressure increases (Fig. 1). The fish communities from intensively exploited sites represent a nonrandom subset (mainly smaller-bodied, lower trophic level taxa) of the species present on unfished reefs. Although we could not directly assess the relative role of fishing and other human effects, we suspect that fishing has been the main proximate cause of the differences between observed and predicted predatory fish biomass, whereas the ultimate causes include coastal development, increased human populations, and economic growth.

The observed total fish biomass in our study varied by nearly 10-fold (55 to 484 g m−2) among the 39 sites; this finding is consistent with other large-scale reef studies across gradients of human impacts. Our range fell within the wider range (15 to 596 g m−2) observed in similar fore-reef habitats of the Caribbean () and elsewhere (note that we did not survey some of the region’s most overfished islands, including Haiti and Jamaica). In the western-central Pacific, for instance, fish biomass gradually increased from 13 g m−2 on reefs of the heavily populated island of Guam to 348 g m−2 on the isolated Kure atoll (), and up to 527 g m−2 on the remote Kingman atoll (). This generalized gradient of fish biomass across large spatial scales is assumed to have been caused largely by spatial variation in fishing intensity because of proximity to human settlements (, ). Because our survey sites were not selected at random (that is, we purposefully surveyed most regions with the least fished reefs), the biomass mean and distribution we documented among sites are not representative of the Caribbean mean/distribution. We suspect that the Caribbean distribution is far more strongly skewed, where a large majority of sites have very low biomass ().

In conclusion, current predatory reef fish assemblages are strongly influenced by both fishing and numerous habitat characteristics. Because there is little spatial coincidence between the natural and anthropogenic factors that influence predator populations, predicting their relative effects at a given location has been challenging. However, our results suggest that not considering the role of natural environmental variation could lead to large underestimates or overestimates of reef-specific carrying capacity for predatory fish assemblages and, hence, the degree of overfishing. We could overlook reefs that provide ideal conditions for predators or greatly overestimate the restoration potential of seemingly predator-depleted sites that may have never supported large populations because of suboptimal physical and biological conditions. Therefore, the carrying capacity for predatory reef fish biomass should be site-specific, and the proposed global unfished biomass derived from remote sites or marine reserves with unique oceanographic features is unlikely to provide an accurate representation of potential carrying capacity conditions in most locations.

MATERIALS AND METHODS

Study sites

Surveys were performed on slope and spur-and-groove fore-reefs with a depth of 10 to 15 m, usually dominated by the corals Montastraea and Orbicella, across 39 sites in The Bahamas, Cuba, Florida (United States), Mexico, and Belize (Fig. 4 and table S2). We selected sites to maximize the range of total fish biomass in each subregion by including reefs inside and outside marine reserves (that is, NTZs, where fishing is prohibited), except at Dry Tortugas, where only a reserve site was surveyed. Four sites in Gardens of the Queen marine reserve in Cuba (Fig. 4 and table S2) were chosen a priori because they were reputed to have relatively high predator biomass with fairly intact fish communities (, ). To minimize seasonal variability, we conducted all surveys during the summer months of May to July 2010–2012 (table S2).

Fig. 4

Distribution of survey locations.

For site abbreviations, survey dates, coordinates, and protection level, refer to table S2. NTZs and minimum fished MPAs are represented with solid symbols. MBR, Mesoamerican Barrier Reef.

Fish abundance

Noninstantaneous underwater visual censuses, with methods modified from Lang et al. (), 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, number, and estimated body size. Fish total length (TL) was estimated and placed into 10-cm size class intervals (), except for individuals with <10-cm TL, for which two 5-cm intervals were used. Precision at size estimation among divers was evaluated before surveys. As transect tape was positioned, a diver swam and counted fish of medium size (5- to 40-cm TL) in (that is, within, passing into, or crossing) a 30 × 2–m belt area, followed by a 15 × 1–m belt to estimate small fish with <5-cm TL (). A second diver counted fish with >40-cm TL in a 50 × 10–m belt to account for more mobile and large-bodied fish (for example, sharks) and minimize overestimation of their densities (). The two smaller transects were contained within the largest transect to create a transect set. Large schools of fish passing through the transect sets were excluded from the data analysis. Each transect set was surveyed in ~15 min, covered the entire visible water column, and was at least 10 m apart from another.

Fish biomass per unit area was the main response variable estimated in our study and used in all the analyses. Biomass was calculated by the allometric length-weight conversion formula W = aTL, where W is body mass in grams and TL is the total length of each fish in centimeters. We used the midpoint of each size class to calculate biomass (), and the parameters a and b are species-specific selected from geographic areas close to our study region (). When these parameters were unavailable, we used estimates for congeneric species of similar morphology and size (table S3). We used fish biomass because it is a comprehensive indicator of fish assemblages status across fishing gradients (, , , ). Fish size is a key functional attribute, and size spectra analyses are used to determine fishing intensity (, ), detect mesopredator release and cascading effects (, ), and even estimate fish abundance in the absence of fishing (). However, fish biomass allowed us to directly compare our results with other studies that estimate unfished biomass (, –, ).

Fish species were assigned to a functional group based on six trophic guilds: apex predators, piscivore-invertivores, invertivores, planktivores, omnivores, and herbivores, following reported dietary information (). For the purpose of this study, we considered “predatory fish” apex predators and/or piscivore-invertivores because both feed mostly on fish. Apex predators were defined as large fish species with the highest trophic level in the community, which consume mostly fish and with very few, if any, predators (). Piscivore-invertivores (that is, mesopredators) feed on fish and invertebrates, invertivores only feed on invertebrates, omnivores consume marine plants and invertebrates, and herbivores only feed on marine plants (table S3).

Covariates

For each reef site, we gathered a preliminary data set of 29 anthropogenic, physical, biotic, and management-related variables known to influence predatory fish abundance (table S4). For a detailed justification of each covariate, see section S4. After an exploratory and discarding analysis, we used the best explanatory variables to predict predatory fish biomass in the absence of humans, that is, no fishing (see “Data analysis”).

Direct and accurate measures of fishing intensity are scarce for our study sites because of the lack of reliable coral reef fisheries data. Therefore, we used human pressure indices, such as human population density and other human-related parameters, as indicators of anthropogenic effects (for example, fishing intensity). Human pressure indices, such as spatial gradients of human population density per reef area or fishers per kilometer of reef, have been used as a proxy for regional fishing intensity in coral reefs, particularly when fishing effort or catch data are not available (, , , , , , , ). For example, the number of people per reef area is strongly correlated with ecological footprints and sustainability of coral reef fisheries across islands of the world (, ). An increasing number of studies have shown that reefs associated with large human population densities have lower abundance and smaller size of fish (, , , , , , ). This approach takes advantage of quasi-pristine reef areas with minimal human disturbance (, , , , ) or of the oldest, well-enforced, and largest fisheries closures as reference sites to evaluate exploitation effects on more disturbed areas (, –). These undisturbed sites should reflect preexploitation levels that can approximate baselines at current conditions (, ).

Anthropogenic variables included coastal development (estimated by electrical power), number of humans, and area of cultivated land (proxy of terrestrial runoff or pollution). These variables were measured within 50 km of each reef site, because this radius of influence has been adequate in detecting human population effects in the region (). Additionally, we measured the number of humans from, and distance to, the nearest population settlement per reef site as a proxy for fish demand and the distance traveled to markets, because both indices have shown to influence local coral reef fisheries (section S4 and table S4) (, ).

Physical and spatial covariates included the following: average and minimum average of sea surface temperature (2002–2011), average oceanic net primary productivity (2002–2012), wave exposure, depth, reef structural complexity, reef area (within 5 and 10 km), distance to deep water (30-m isobath), distance to reef breaks, and distance to mangrove (section S4 and table S4). Biotic factors included mangrove perimeter (within 5 and 10 km), coral cover, macroalgal cover, gorgonian abundance, and biomass of lower trophic fish groups. Reef area and mangrove perimeter were calculated at multiple scales to determine the influence of landscape extent on predatory fish (section S4 and table S4). Management-related variables included protection level for each site (none, MPAs, and NTZs), reserve size and age, and poaching levels (low or high) inside reserves (section S4 and table S4). For detailed descriptions and measurements of each covariate, refer to section S4.

Data analysis

To explore among site variation of predatory fish biomass in relation to strict protection (that is, reserve and nonreserve) and country, we used a linear mixed-effects model in which fish biomass was explained by those two factors and grouped by sites. We analyzed differences between factors using Bonferroni-corrected pairwise comparisons.

To select the covariates that better explain the variation in predatory fish biomass, we first evaluated their collinearity using a Spearman’s rank (rs) correlation matrix for all sites and for sites within marine reserves (table S5). Several covariates were sufficiently correlated (−0.5 > rs > +0.5) to compromise interpretation when modeled together (). For example, reef area (rs = 0.83) and mangrove perimeter (rs = 0.93) were highly correlated within 5 and 10 km, as were the log values among most of the human-related variables (rs > 0.5) (see table S5 for other correlations). Thus, we first ran single-variable generalized linear models with related covariates (for example, human-related) to examine the best supporting covariates using the weights of AICc (table S6). Improvement in fit was evaluated with analysis of deviance among models ().

We created two sets of global models: set A, for all reef sites considering protection level, and set B, for the subset of sites within marine reserves (Table 1). This approach allowed us to assess the effect of protection level among sites (model set A) and analyze the variation of predatory fish biomass within reserves, where fishing is theoretically absent but enforcement varies among sites (model set B). Within each model set, we modeled separately those covariates that were correlated (table S6), eliminating factors that did not improve model fit. To verify the lack of multicollinearity among covariates, we calculated the variance inflation factors (VIFs) after fitting the models. We sequentially removed and modeled separately each covariate for which the VIF value was above 2 ().

All the human-related variables, except “distance to population centers,” explained some of the variation of predatory fish biomass among sites in the single-variable models (table S6). The “log of coastal development within 50 km” (hereafter coastal development) yielded better goodness of fit (that is, the lowest AICc and the highest weights) when considering all sites. However, for the subset of sites within marine reserves, the “log of humans within 50 km” (hereafter human population) showed higher weights for total predators and piscivore-invertivore biomass, whereas the “log of cultivated land within 50 km” (hereafter cultivated land) had the highest weights for apex predators. We selected “minimum monthly sea surface temperature” (hereafter temperature), “reef area within 5 km” (hereafter “reef area”), and “mangrove perimeter within 5 km” (hereafter “mangrove”) because they had the highest AICc weights (table S6). We discarded physical variables such as “wave exposure,” “depth,” “distance to deeper water,” “minimum distance to mangrove,” “distance to reef breaks,” and the biotic variable “macroalgal cover” because they did not contribute to model fit in exploratory models. Different combinations of noncorrelated variables were considered candidate predictors for predatory fish biomass.

We evaluated the relationship between the selected set of variables and the biomass of predatory reef fish (apex predators, piscivore-invertivores, and total predators) with GLMMs () fitted by maximum likelihood (Laplace approximation), in which reef sites, region, and year of survey were nested and coded as random effects and the explanatory variables as fixed effects. Exploratory analysis showed that predatory fish biomass did not change between years of survey. Thus, this random structure was used to nest the variance and account for the spatial autocorrelation observed in the raw data (, ). The biomass of total predators and piscivore-invertivores was modeled with a Gaussian distribution and log link. For apex predators, we used a zero-inflated model with Gaussian error structure and log link to account for the excess of zeroes in the biomass distribution. The biomass of all fish groups, including response and predictors, was log10(x + 1)–transformed to improve homogeneity of variance and model fit. Numerical covariates were standardized and centered (mean of 0 and SD of 1) to aid in model comparisons. Meaningful interactions and quadratic terms were included in exploratory models.

A multimodel inference approach and model averaging based on AICc weights (ΔAICc < 2, where ΣAICc weights > 0.95) were used to select the variables included in the best models within each set (Table 1) (). For the final models, we calculated a pseudocoefficient of determination for GLMMs (pseudo-R2), interpreted as the variance explained by both fixed and random factors, that is, the entire model (). Homogeneous and normal distribution errors of final models were confirmed in the plot of residuals against fitted values and by using the normal scores of standardized residuals deviance, respectively (fig. S6) (). Spline spatial correlograms were plotted to corroborate that the final model residuals were independent and not spatially autocorrelated (fig. S7) (). All analyses were performed in R v.3.03 () using the packages lme4 v.0.99-2 () and glmmADMB v.0.8.0 () for GLMMs and MuMIn v.1.15.6 () for model averaging.

Carrying capacity prediction

Finally, using the best explanatory models for total predatory fish, we predicted the expected biomass range at each reef site in the absence of humans by setting human-related variables to zero and categorizing all sites as NTZs (that is, no fishing) (section S5). This approach was used before to estimate a potential range of unfished sharks and reef fish biomass across islands of the central and western Pacific (, ). Here, the predicted values are potential biomass of predatory fish at each site under current conditions. We assumed that the trend between human impacts and predatory fish biomass can be linearly interpolated from relatively low levels of human effects back to zero. This may underestimate potential carrying capacity because predators are vulnerable to even low levels of fishing in coral reefs (, ). In addition, we assumed that interaction between predators and lower trophic groups responds linearly and would extrapolate as such without humans. Finally, we assumed no trophic cascades of predators on prey within fisheries closures (). We back-transformed the values of predatory fish biomass and presented the quantile median values in a summary figure. Predicted means and 95% confidence intervals per site were obtained by bootstrapping (section S5).