Non-reef habitats such as mangroves, seagrass, and macroalgal beds are important for foraging, spawning, and as nursery habitat for some coral reef fishes. The spatial configuration of non-reef habitats adjacent to coral reefs can therefore have a substantial influence on the distribution and composition of reef fish. We investigate how different habitats in a tropical seascape in the Philippines influence the presence, density, and biomass of coral reef fishes to understand the relative importance of different habitats across various spatial scales. A detailed seascape map generated from satellite imagery was combined with field surveys of fish and benthic habitat on coral reefs. We then compared the relative importance of local reef (within coral reef) and adjacent habitat (habitats in the surrounding seascape) variables for coral reef fishes. Overall, adjacent habitat variables were as important as local reef variables in explaining reef fish density and biomass, despite being fewer in number in final models. For adult and juvenile wrasses (Labridae), and juveniles of some parrotfish taxa (Chlorurus), adjacent habitat was more important in explaining fish density and biomass. Notably, wrasses were positively influenced by the amount of sand and macroalgae in the adjacent seascape. Adjacent habitat metrics with the highest relative importance were sand (positive), macroalgae (positive) and mangrove habitats (negative), and fish responses to these metrics were consistent across fish groups evaluated. The 500-m spatial scale was selected most often in models for seascape variables. Local coral reef variables with the greatest importance were percent cover of live coral (positive), sand (negative), and macroalgae (mixed). Incorporating spatial metrics that describe the surrounding seascape will capture more holistic patterns of fish-habitat relationships on reefs. This is important in regions where protection of reef fish habitat is an integral part of fisheries management but where protection of non-reef habitats is often overlooked.

Study Location

This study was conducted around Siquijor Island in the Visayan region of the Philippines (Fig. 1a). Shallow water benthic habitats of Siquijor include macroalgal beds, mangroves, and seagrass beds of varying spatial extent adjacent to fringing coral reefs. 

Fish and Habitat Surveys

Surveys of reef fish and benthos were conducted in April – July 2016 at eight locations around Siquijor Island (Fig. 1a), with paired NTMR and control (open to fishing) sites, totalling to 16 sites. Along a 50-m by 5-m transect, large mobile reef fish (>10 cm TL) were counted and sized to the nearest centimetre. On the return swim, smaller (≤ 10 cm TL) reef fish species were recorded within a 2-m width. Biomass of fishes was calculated using published length-weight relationships (Kulbicki et al 2005). For benthic surveys, substratum was identified at 50-cm intervals along the 50-m transect and was classified based on substrate (rock, sand, rubble, coarse sand) and benthic cover (abiotic, crustose coralline algae, epilithic algal matrix, macroalgae, soft coral, hard coral, other) (Table 1). Macroalgae and soft coral were identified to genus when possible. Hard coral was identified to genus and classified into growth form (fragile, robust). The ‘other’ category included sessile invertebrates such as sponges, tunicates, and gorgonians. Structural complexity was estimated visually on a 0-5 scale following methods used in Wilson et al. (2007). 

Habitat Mapping

 Images from the GeoEye and PlanetScope satellite sensors were acquired from the Digital Globe Foundation, and Planet, respectively. The GeoEye satellite provides a spatial resolution of 1.84-m and Planet provides a 3-m resolution, both across four spectral bands of blue, green, red, and near-infrared (NIR) (Fig. 1b). Both sensors were necessary to acquire complete coverage of the island. Pre-processing of imagery was conducted using the software ENVI (v. 5.3, Harris Geospatial Inc.). Band ratios were calculated to provide additional unique spectral signatures for benthic habitat classes (Phinn et al. 2012, Roelfsema et al. 2013). Band ratios were: blue to red (B/R), blue to green (B/G), and red to NIR (R/NIR). After pre-processing, classification of imagery into habitat types was conducted using the maximum likelihood classification tool in ArcGIS, v. 10.4.1. Feature classes were a combination of biotic and geomorphological features: seagrass meadows, macroalgal beds, reef flat, reef crest, reef slope, lagoon, sand, mangrove forest, and beach (Fig. 1c). Georeferenced habitat data points (n=500) collected in-situ in 2016-2018 informed the maximum likelihood classification, with 70% of points used for training, and the remaining 30% used for validation of the classified map. The map was then manually reviewed and edited for obvious errors, smoothed using the majority filter in ArcGIS, and converted to polygons for spatial analysis. Map validation identified 72% accuracy of habitat classification using the maximum likelihood method.

Spatial Analysis

Fish and benthic survey locations were overlaid onto the classified habitat map to calculate spatial statistics of the seascape surrounding each site (n = 16). Adjacent habitats used for spatial analysis were seagrass, macroalgae, sand, reef flat, and mangroves. For each location, distance to the nearest habitat type was measured using edge-to-edge distance between survey sites and each habitat. Buffer zones surrounding each survey site were calculated at three different spatial scales (250, 500, 1000-m) (Fig. 1c). Buffers were clipped by shore and deep-water features to only represent shallow water habitat. The proportion of each habitat within each buffer zone was calculated as the area of habitat divided by the total area of the clipped buffer. These data were then incorporated with the benthic survey data on coral reefs for further analysis (Table 1).  Global Moran’s I was calculated for the 500-m habitat spatial scale to evaluate any potential spatial autocorrelation. Spatial data was not significantly spatially autocorrelated for the 500-m scale (Moran’s I = 0.370, p = 0.24).

Statistical Analysis

Boosted regression trees (BRT; Elith et al. 2008) were used to evaluate how benthic habitats at different spatial scales affected coral reef fishes using the gradient BRT method from the gbm package. Fish groups were analysed in terms of density and biomass, or presence/absence, using Poisson, Gaussian, and Bernoulli distributions, respectively. Presence/absence was used for species groups with too few observations for density and biomass analysis (Lutjanidae and Serranidae). In total, 32 BRT models were run on fish groups with the greatest number of observations at the family level: Labridae (wrasses, excluding parrotfishes), Lutjanidae (snappers), Serranidae (groupers), Pomacentridae (damselfishes), Chaetodontidae (butterflyfishes), and Acanthuridae (surgeonfishes) (Table 2). Parrotfishes (Labridae, subfamily Scarinae) were run at the level of genus for two different feeding-type groups, Scarus and Chlorurus, where Scarus are scrapers and Chlorurus excavators. Hipposcarus was included in the ‘Scarus’ group and Cetoscarus was included in the ‘Chlorurus’ group based on their feeding modes. Models for juvenile reef fish density were only possible for wrasses, and the parrotfish groups Scarus and Chlorurus, due to the lack of juveniles observed from other families. Fish groups were also separated by coral reef zones, i.e. reef crest and slope.

To identify the scale at which reef fish responded to the seascape, a BRT was run for each adjacent habitat type at all three spatial scales (250, 500, 1000 m) for each response variable. The ‘best’ scale for each habitat type was selected as the radius with the highest relative importance, and only that scale was included for further analysis. Variables with correlation values greater than 0.8 (e.g. hard coral, fragile coral, robust coral) were run in a BRT, and only the variable with highest relative importance was selected for the remaining analysis. Full models were then run with these pre-selected variables with an interaction depth of 3 and bag fraction of 0.75 using the gbm.step method in the gbm package, and were calibrated for best results by altering the learning rate to achieve the optimal number of iterations between 1000-10000 trees, based on a 10-fold cross-validation procedure. The gbm.simplify process was used to reduce the number of variables by an iterative backwards stepwise removal of the least influential variables using k-fold cross validation until the change in predictive deviance was minimized. The simplify process selected the nine most influential variables, and NTMR status was the tenth variable to evaluate any reserve effect. To account for stochasticity and incorporate uncertainty values for relative importance, models were bootstrapped (sampling with replacement) 100 times. Error in relative importance and deviance explained values were measured by 95% confidence intervals from the bootstrapping process. Cross validation deviance (CV deviance) was calculated by subtracting the CV deviance from the null deviance and dividing by the null deviance. Mean relative importance was used as an indicator for variable importance. Because models had 10 variables, relative importance values greater than 10% were considered influential as they were selected more frequently than expected by chance. The mean relative importance was summarised only for influential variables (>10% relative importance) and compared between variables categories (local reef vs. adjacent habitat) (Table 1). Here, we define “local reef” as the small-scale benthic habitat characteristics of a coral reef, whereas “adjacent habitat” describes larger scale spatial metrics of multiple habitat types across a seascape. Wilcoxon ranked tests for non-parametric data were used to compare the mean relative importance between local reef and adjacent habitat categories across all models, at the level of reef zones (crest and slope), fish life stages (juvenile and adult), and for each fish group.

Transect_Summ_data_for_Dryad.xlsx:
Data for analysis, with descriptions of variables.

Bootstrap_BRT_full_dryad.Rmd:

.Rmd file detailing variable selection, gbm process, and boostrap example when running data for one respone variable.

BRT_BOOTSTRAP_final.Rmd: 

.Rmd file for bootstrap to run through multiple models at once. Loops through bootstrap and also loops through mulitple models.