Open Collections

Description

Abstract

Coral reef ecosystems are widely threatened by global change, yet the cumulative impacts of multiple interacting stressors remain difficult to quantify over space and time. We evaluate how long-term artisanal fishing effort, blast fishing, human population density, and marine protected areas (MPAs) interact with biophysical and seascape variables to influence the spatial distribution of corals in a biodiverse, but heavily impacted ecosystem.

To address these challenges, we combined satellite habitat mapping, indigenous and local knowledge, and generalized linear mixed-effects models to assess how stressors and seascape characteristics shape the occurrence of living coral. This integrative approach allowed us to capture both ecological and social processes at an ecosystem scale.

Coral was the dominant habitat in only 30% of the study area. The strongest predictor of coral distribution was seascape configuration: corals were more likely to be in compact reef patches. Coral presence was also positively associated with MPAs (18% higher probability inside MPAs) and depth (8% increase in probability of corals per 2.25 m). Increasing distance to seagrass was associated with higher coral probability (6% per 100 m), but this effect diminished at greater distances, reflecting a nonlinear quadratic response.

Disturbance effects were complex. Coral probability declined 6% for each 35 days of blast fishing (fishing with explosives), with an antagonistic interaction with human population density: impacts were strongest in low-population areas (−7.6%) and weaker in high-population areas (−2.4%). Total fishing effort (excluding blast fishing) reduced coral probability by 3.2% for every ~530 fishing days at a site, with a notable 10-year lag, highlighting delayed ecosystem responses.

Our findings emphasize the importance of long time series and the need to account for lagged effects of fishing, which may otherwise be underestimated. They also underscore that conservation outcomes depend on both managing harmful stressors and enhancing beneficial seascape features. Key priorities include reducing spatial overlap of destructive activities, protecting reef configurations that support coral persistence, and addressing stressors with delayed impacts. The approach developed here—integrating spatial, ecological, and social data—offers a framework for adaptively managing coral reef ecosystems that are changing under interacting pressures.

Study site

The Philippines, located in the global center of marine biodiversity, supports 21,983 km2 of coral reefs, and is a global priority for conservation due to the high threats to the region (Halpern et al., 2008; Roberts et al., 2002) and a relatively high resilience to coral bleaching (Sully et al., 2019). The Danajon Bank (10˚15’0’N, 124˚8’0’E) is a double barrier reefs in the central Philippines that sits off the northern edge of Bohol – a province characterized by small-scale farming and high levels of erosion (Provincial Government of Bohol, 2011). Several characteristics of the Danajon Bank potentially enable it to resist anthropogenic stressors. Such characteristics include high biodiversity and low incidence of coral bleaching (excluding the 1998 bleaching event) (Bayley et al., 2019; Philippine Coral Bleaching Watch, 2013; Sully et al., 2019), although data for many stressors are sparse or at coarse scales (Table S1; Data Sources). Located at the southern end of the sheltered Camotes Sea, the Danajon Bank’s oceanography is characterized by a well-mixed water column, alongshore flow, weak tidal circulation, and the absence of large waves (Villanoy et al., 2006). 

The Danajon Bank lies adjacent to Cebu City (the second largest metropolitan area in the Philippines, 3.1 million people), but is itself located in a rural region, which struggles with extreme poverty (Provincial Government of Bohol, 2011). Small-scale fisheries are based on the mainland of Bohol and on small cays, and have exclusive fishing rights to coastal areas (Christie et al., 2006). These fishers use many fishing gears and target diverse species of marine life (Kleiber et al., 2014; Selgrath et al., 2017). Since 1998, small-scale fisheries have been co-managed and many communities have established locally-enforced no-take marine protected areas (MPAs), including 12 MPAs in our study area (Pajaro, 2010; Selgrath et al., 2018). MPAs in the region were established by local fishing communities and NGOs, and frequently placed in areas that had been degraded by fishing prior to their protection (Hansen et al., 2011).  Despite these management efforts, heavy fishing pressure, destructive fishing, and high population densities have depleted marine life and degraded much of the marine environment (Bayley et al., 2020; Magdaong et al., 2014; Selgrath et al., 2016).

Overview

This spatial analysis focused on a 19km by 22km area in the central Danajon Bank (hereafter ‘Study Area’). In the Study Area, we evaluated the relationship between reef state (presence of living coral vs rubble) and 25 spatially explicit variables (Table S1; Data Sources). Variables described anthropogenic and biophysical attributes of the reef system. These social-ecological attributes were derived from remote sensing imagery, participatory mapping with local fishers, and publicly available data sources (Table S1; Data Sources). Using mixed, logistic models, we considered additive, antagonistic, and synergistic effects of stressors to understand their influence on the system. 

For this analysis, we used a binary measurement of reef state, classified as (a) ‘coral-dominated’, or (b) ‘rubble dominated’ (see Marine habitats below for classification details, hereafter ‘coral’ or ‘rubble’). Using a binary variable as a proxy for reef functioning represents the inherent trade-offs between highly detailed surveys of small areas and the coarser mapping of whole ecosystems. 

Spatial datasets Marine habitats

We created maps of marine habitats of the Study Area by classifying high spatial resolution (2 m) multi-spectral satellite images (WorldView2, dates: 2010/05/10 and 2012/20/04) using object-based image analysis (OBIA). These two dates were the earliest times with clear WorldView2 images of the Study Area, as the satellite was launched in October 2009. Mapping was restricted to shallow areas (approximately 15 m depth). Geomorphic zones (e.g. reef crest, reef flat) were identified and further classified into 16 benthic (seafloor) habitat types or classes (Table S2). 

For this analysis, we simplified coastal habitats into 5 classes: (living) coral, rubble, sand, seagrass, and mangroves (Table S2). The rubble-dominated category included degraded rubble, deal coral, and macroalgae. Since the Danajon Bank contains areas with mixed habitats, the rubble-dominated class included areas containing rubble mixed with other habitats (e.g., sand, seagrass). During map classification, we used an approximately 30% cover threshold from field survey data to classify reefs as ‘living coral’ during the classification process. As the satellite data map did not capture habitats beds in turbid shallow waters or in deeper waters, we supplemented the habitat maps with ILK information about habitats that we documented during participatory mapping interviews with fishers. See Selgrath et al., 2016 for details on habitat mapping. Using independent data from a long-term monitoring program (Project Seahorse Foundation/ZSL Philippines) we created a confusion matrix to evaluate map accuracy and updated three categories to improve classification accuracy (Table S3). 

Fishing and other anthropogenic pressures

Fishing effort. Fishing is one of the greatest threats to coral reef ecosystems (Burke et al., 2011), but longitudinal spatial measures of fishing effort on reefs are exceedingly rare (Selgrath, Gergel, et al., 2017). We developed maps of fishing effort and fishing gear use during semi-structured, participatory mapping interviews. We interviewed 391 randomly sampled male fishers from 23 randomly selected villages (approximately 7% of fishers from 50% of villages). Villages were located both within the Focal Area, and up to 10km outside of the Focal Area. Of the total respondents, 295 individuals fished inside of the Focal Area. Before conducting interviews, we obtained written consent from local mayors and/or councils, and oral consent from village officials and fishers. Research methods, including consent procedures, were approved by the University of British Columbia’s Human Behavioural Research Ethics Board (H07-00577). In the mapped area, fishers reported using 88 fishing gears (e.g., hand-line, bottom-set gillnets, skin diving, blast fishing). We estimated that during 2010 approximately 8,000 men fished full- or part-time inside the 418 km2 study area (Selgrath et al., 2017).

To evaluate fishing effort, we digitized fishing grounds from interviews using heads up digitizing in ArcGIS 10.1 (Environmental Systems Research Institute, Redlands, California). Maps depicted fishing effort in 20 m by 20 m grid cells. For each year, we integrated individual respondent’s fishing maps to map cumulative spatial fishing effort (SpEFti) as the estimated number of fishing days (E) by all fishers from study villages (F) in year (t) in grid cell (i). Estimating total effort allowed us to account for changes in the sample size of respondents over time. We based the estimated number of fishers (F) in year (t) on estimates of demographics changes in study villages. We also conducted an area accumulation analyses to ensure that we interviewed a representative number of fishers (see Selgrath & Gergel, 2019 for the methods and detailed results of the area accumulation curve analysis).

To evaluate the long-term influence of fishing effort, we normalized pixel values of the yearly fishing effort maps from 0-1 relative to the maximum effort from all years (maximum fishing days per year in grid cell (i) = 5,546 days) (Halpern et al., 2015) where normalized SpEFti is the spatial fishing effort (SpE) for all fishers from study communities (F) in year (t) in grid cell (i) scaled to the maximum fishing effort by all fishers (F) from all cells (I) in all years (T)

*           normalized SpEFti = SpFti/MaxSpEFTI                                                         (1)

We considered the influence of normalized fishing effort (normalized SpEFti) over three time periods:

(1) Contemporary effort (2010): cumulative fishing days (SpE) by all fishers from study communities (F) in year (2010) in grid cell (i)

*           Contemporary Effort2010= SpEF2010i                                                        (2)

(2) Cumulative effort (1980-2010): cumulative number of fishing days (SpE) by all fishers from study communities (F) in two to four years (t1, …, t2010) in grid cell (i)

*           Cumulative Effortt1:2010= SpEFt1i+…+SpEF2010i                                    (3)

(3) Lag effort (1980-2000): cumulative number of fishing days (SpE) by all fishers from study communities (F) in three years (1980, 1990, 2000) in grid cell (i)

 *           Legacy Lag Effort1980:2000= SpEF1980i+SpEF1990i+SpEF2000i             (4) 

We created fishing pressure maps using the same methods for blast fishing using explosives, which we hypothesized would have a large impact on the probability of an area supporting living coral (see Do multiple stressors affecting corals have additive, synergistic, or antagonistic impacts? below).  Blast fishing effort was normalized relative to the maximum blast effort from all years (maximum blast fishing days per year in grid cell (i) = 336 days).  

MPAs, human population, and markets. We assessed the influence of protection by considering sites located inside and outside of existing MPAs. We used population information collected from village census data to evaluate the influence of four anthropogenic variables: population density; distance to the regional fish market in Cebu City; distance to towns; and distance to fishing communities (Table S1). Several ocean stressors (e.g. nutrient loading, trampling) are correlated with the population densities of adjacent communities (Mora, 2008; Rangel-Buitrago et al., 2024). We therefore used distance decay with a square root decay to calculate the risk from population density that decreased with distance to the communities (i.e., an environmental risk surface; Table S1) (McPherson et al., 2008).

Biophysical attributes

Habitat configuration. From the habitat map described above, we derived landscape pattern indices describing the arrangement of habitats in the seascape (Table S1)(Wedding et al., 2011). We emphasized class-level indices describing coral and rubble patches as well as their relationships to other habitat classes (Table S1; Fig. 2)(McGarigal & Cushman, 2002). For points within coral and rubble patches, metrics included: nearest neighbor distance between points and neighboring patches; patch area; patch edge length; and patch compactness (incorporating edge:area ratio and patch size - inverse of patch shape index described in McGarigal et al., 2012). Because reefs exhibit zonation patterns (Hoegh-Guldberg et al., 2011), we also determined the distance of points within coral/rubble patches to mangroves and seagrass habitats (i.e., the distance between the coral/rubble patch and the edge of the nearest mangroves and nearest seagrass beds. We divided distance by 100 m to look at changes per 100 m of distance (Gelman & Hill, 2007).

Bathymetry, rivers, zonation, aragonite saturation, thermal stress. We obtained information for biophysical attributes using existing data sources (Table S1; Fig. 2). First, we created a bathymetry layer by interpolating depth soundings from a NAMRIA nautical chart. Spline interpolation creates a continuous raster surface from sampled point values using a two-dimensional minimum curvature spline technique. Second, we identified river mouths using WorldView-2 satellite images and characterized the seascape based on distance to rivers. Third, we classified the study area into six ecological zones modified from Hansen et al., 2011: (1) inshore coastal zone of turbid waters, extensive mangroves, and some larger communities supporting rice farming and mining; (2) terrestrial islands that support extensive mangroves and farming; (3) coastal inner reef zone near terrestrial islands with moderately clear waters and mangroves; (4) an inner reef zone with intermediate characteristics of the coastal and outer reefs; (5) an offshore, outer reef zone with clear waters, extensive seagrass beds, scattered cays adjacent to the channel between the Danajon Bank and Olango Island; and (6) an offshore, outer reef zone with clear waters, extensive seagrass beds that is adjacent to the Camotes Sea. Fourth, we used data from regional models of aragonite saturation and past thermal stress (1998-2007) (Burke et al., 2011). Calculations used R (R Core Team, 2025) and ArcPro 3.5 (Environmental Systems Research Institute, Redlands, USA).