4.3.1. Site characteristics and sampling design
A total of 182 reefs across 63 locations within 17 different Pacific Island countries and territories were surveyed once between 2003-2008 (Figure 4.1), within the framework of the Pacific Regional Oceanic and Coastal Fisheries Development Programme (PROCFish/C/CoFish) under the auspices of the Pacific Community (SPC) (for sampling methodology and rationale, see Appendix Bi and Pinca et al., 2009). Reef fish and benthic communities were monitored at reefs within each site, and landing surveys quantified fisheries catch. Fishing grounds for each site were delineated from information given by local fishers. Total reef area (km2) within each site’s fishing ground was then quantified from satellite images, allowing for subsequent calculation of human population (total population within villages with access to the fishing ground) and annual fisheries catch relative to reef area (i.e. local human density and relative finfish catch).
4.3.2. Field surveys
Reef fish communities were measured using the distance-sampling underwater visual census method along 50m transects (described in Labrosse et al., 2002). Data on abundance and body length were recorded to species level for herbivorous fish. Counts were converted to biomass (g m-2) for each species from established length-weight relationships (Kulbicki et al., 2005) and averaged across transects within each reef habitat type. Benthic cover data was obtained using the medium-scale approach (described by Clua et al., 2006). The method is based on a semi-quantitative description of 10 quadrats
LOCAL HUMAN IMPACTS AND CORAL REEF BENTHIC ASSEMBLAGES
68
of 25 m2 (5 x 5 m) laid down on each side of the 50 m transects used for reef fish surveys (500 m2 per transect in total). Surveyors first recorded abiotic and live coral substrates (e.g.
sand, rubble, rocky slab and boulders, and hard coral - live, bleached and long dead, with live coral divided into various morphologies e.g. branching, massive). Each component was quickly estimated using a semi-quantitative scale, estimating in units of 5% and ranging from 0 to 100%. Secondly, groups growing on the substrate (e.g. fleshy algae, cyanobacteria, turf algae, crustose coralline algae/CCA) were recorded using the same semi-quantitative scale. Joint fish-benthic transect replication varied among reefs (n = 3 to 47). Transect data were pooled for each reef habitat type within each location.
Figure 4.1. Map of 63 survey sites; temporal structure shown by colours related to survey year.
4.3.3. Intrinsic reef attributes
Dominant benthic groups (mean composition >10% of benthos) were dead coral (incorporating rubble, boulders and pavement), live hard coral, turf algae, fleshy algae and crustose coralline algae (CCA). Pairwise relationship tests (corvif function in R –Zuur et al., 2009) between groups established no collinearity (R2 < 0.5). Among-reef variation of benthic assemblages was then explored using a principal components analysis (PCA) based on Euclidean distances (built-in prcomp function in R - R Development Core Team 2013). A hierarchical clustering of the dominant benthic groups was produced with the same Euclidean distance matrix as for the PCA using pvclust package in R (Suzuki &
Shimodaira, 2015). Significance values were calculated by 10,000-fold multi-scale bootstrap resampling, with independent clusters considered for significance values greater than 0.95. The first (PC1) and second (PC2) principal components of the PCA were then extracted for each site to provide one-dimensional multivariate representations of the benthic communities. Selected response parameters included arcsine transformed percentage cover of live hard coral, turf algae and fleshy algae, as well as the extracted PC1 and PC2 (Table 4.1). To explore patterns in morphological composition of the coral community, the ratio of massive to branching morphologies (i.e. higher values indicate relatively more massive corals) was used as an additional benthic response parameter.
69
Table 4.1. Benthic response variables and ranges throughout all reefs.
4.3.4. Environmental predictors implicit for reefs
Model predictors (Table 4.2) included a variety of factors that were either collected during PROCFish surveys or extracted from public repositories. Depth was averaged over transects within each habitat at each site, and habitat represented reef habitat type (e.g.
back, coastal, lagoonal, outer; Appendix Bi). Collinearity between depth and habitat precluded these predictors being used in parallel, thus depth was incorporated into models focused on outer reefs and habitat in those including different reef habitats. Due to natural latitudinal differences in reefs (Hughes et al., 1999; Harriott & Banks, 2002), latitude was represented by degree distance from the equator without differentiating between north and south (0 to 23.9°). Average values for degree heating weeks (DHW - sum of previous 12 weeks that thermal stress anomalies >=1°C) at each site (n = 55; no data for some locations) were extracted from the NOAA Coral Reef Thermal Anomaly Database (CoRTAD version 4 – Casey et al., 2012), at a spatial resolution of 4 km. Thermal stress anomalies refer to weekly sea surface temperature minus maximum weekly climatological sea surface temperature (Casey et al., 2012). Data were averaged over the 12 years prior to each specific site’s survey date, allowing all sites to encompass the strong
Response Description Range
CCA Average benthic cover (%) of crustose coralline algae
0.0 to 43.5 dead coral Average benthic cover (%) of dead coral;
including rubble, boulders and pavement
5.5 to 61.2 fleshy algae Average benthic cover (%) of fleshy algae 0.0 to 51.8 turf algae Average benthic cover (%) of filamentous algal
turfs 0.0 to 45.6
live hard coral Average benthic cover (%) of live hard coral 6.1 to 65.1 pc1 Extracted first principal component from
principal components analysis of key benthic groups; explaining a trend from fleshy algae and dead coral dominated reefs (negative values) to live hard coral and CCA (positive values)
-3.4 to 3.2
pc2 Extracted second principal component from principal components analysis of key benthic groups; explaining a trend from low (negative values) to high (positive values) turf algal coverage
-2.9 to 2.4
coral morphologies Proportion of branching to massive coral morphologies within the live hard coral community; increasing values represent relatively more massive morphologies
0.0 to 70.2
LOCAL HUMAN IMPACTS AND CORAL REEF BENTHIC ASSEMBLAGES
70
1998 El Niño event as well as the 2000 and 2001 bleaching events in the region (Appendix Bi), and also considering that reefs have been found to recover from acute disturbances within this timeframe (Gilmour et al., 2013). Strong negative collinearity between latitude and DHW (R2 > -0.9, Appendix Bii) precluded including both predictors in combination.
Because of the complete data for latitude, this predictor was selected to represent both factors. Storm exposure was extracted from the NOAA IBTrACS-WMO data (Knapp et al., 2010a, 2010b) within ArcMAP 10.4 (ESRI, 2011). After projecting sites using the Behrmann projection, the number of storms (category one through five) passing within a 50 km buffer of each site was extracted. Based on the same justification as for DHW, storm exposure data were also confined to 12 years prior to the survey date for each individual location.
From socioeconomic surveys, the predictor relative local human density was selected since it exhibits positive collinearity with both relative finfish catch and number of boats relative to reef area (R2 > 0.65), and was thus assumed to represent relative fishing pressure as well as associated effects of human populations on coastal water quality. All herbivorous fish species encountered during visual surveys were classified into functional groups (Appendix Biii) according to Green & Bellwood (2009), enabling biomass (g m-2) of grazers/detritivores, browsers, scrapers/small excavators, and large excavators/
bioeroders to be incorporated as model predictors.
Due to large variation in the scales of different terms, predictors underwent z-transformation to allow appropriate comparisons between the effect sizes (Zuur et al., 2009). Prior to modelling, pairs plots were assessed for collinearity between model terms (corvif function – Zuur et al., 2009; Appendix Bii). Multi-collinearity was additionally checked using the generalised variance inflation factor (GVIF) function in R (car package – Fox & Weisberg, 2011) where values >3 would suggest collinearity and were not observed.
4.3.5. Defining level of human impact
We categorised reefs into those exposed to low and high local human impact according to relative local human density. Low (n = 29 outer reefs, n = 99 reefs encompassing all habitat types) and high (n = 33 outer reefs, n = 83 reefs encompassing all habitat types) local human impact was established around a threshold of 25 people km-2 reef. This threshold was previously identified as a breaking point within the same data-set after which phylogenetic diversity of parrotfishes was significantly reduced (D’Agata et al., 2014). Relative local human density was also collinear (R2 = 0.7) with relative finfish catch (tonnes reef fish km-2 reef year-2). As an additional analysis allowing us to explore more specifically the impact of fisheries pressure, we categorised reefs into those exposed to low (n = 35 outer reefs, n = 107 reefs encompassing all habitat types) and high (n = 27 outer reefs, n = 75 reefs encompassing all habitat types) fishing pressure using a threshold of 5 tonnes reef fish km-2 reef year-1. This value reflects a conservative estimation of maximum sustainable yield (MSY) of coral reef fisheries (Munro, 1984; Dalzell & Adams, 1996; Newton et al., 2007). These two separate classifications of local impacts tested a slightly different set of reefs focused on either human densities or fishing pressure. However, because of similarities between the two factors (R2 = 0.8), and difficulties determining whether low catch indicates low fishing pressure or overexploited resources, results according to relative finfish catch are confined to the appendix (Appendix Biv). All references hereafter to local human impact refer to those classified according to relative local human density.
71
Table 4.2. Model predictors and ranges among all reefs.
Predictor Description Range
browser biomass Biomass (g m-2) of browsers 0.1 to 73.2 degree heating
weeks (DHW)
Measure of cumulative thermal stress – sum of previous 12 weeks where thermal stress
anomaly ≥1°C; value averaged over 12 years preceeding survey; negatively collinear (R2 = 0.9) with latitude; only available for n=55 sites
0.6 to 3.5
depth Average depth (m) of transects 1.25 to 10.5
excavator biomass Biomass (g m-2) of large excavators/bioeroders 0.0 to 1103.4 grazer biomass Biomass (g m-2) of grazers/detritivores 1.1 to 161.0 habitat Reef habitat type; grouped into four
geomorphological structures: sheltered coastal reef, intermediate lagoon reef (patch reef inside lagoon), back reef (inner side of outer reef) and outer reef (exposed reef) latitude Degrees (°) distance from equator, not
differentiating between north and south.
Negatively collinear (R2 = 0.9) with degree heating weeks (DHW)
0.0 to 23.9
local human density
Number of people km-2 reef, positively collinear (R2 = 0.8) with relative finfish catch
1.3 to 1705.0 relative finfish
catch
Annual reef finfish catch (tonnes) km-2 reef year-1; positively collinear (R2 = 0.8) with local human density
0.1 to 78.2
scraper biomass Biomass (g m-2) of scrapers/small excavators 2.5 to 151.3 storm exposure Total number of storms (category 1 to 5)
passing within 50 km of site within previous 12 years
0 to 14
4.3.6. Model design and selection
To test which predictors show the strongest association with the various benthic response variables, a series of models were constructed (Table 4.3). Two predefined models involved (i) spatial and physical predictors (i.e. latitude, storm exposure, reef depth), and (ii) local feature predictors available from surveying each site (i.e. biomass of herbivorous fish functional groups, local human density, reef depth). Reef depth was included in both models to investigate the respective effects when accounting for this intrinsic factor which
LOCAL HUMAN IMPACTS AND CORAL REEF BENTHIC ASSEMBLAGES
72
can mediate the impacts of both kinds of predictors, for example storms, temperature anomalies and fishing pressure (Tyler et al., 2009; Bridge et al., 2013). Finally, (iii) best-fit1 models were determined from model selection techniques using the MuMIn package in R (Barton, 2016) based on AIC values corrected for finite sample sizes (AICc). From a model containing all predictors (i.e. Table 4.2) the dredge function runs all possible combinations and returns various models ranked from best to worst according to AIC. The function also returns a value between zero and one for each predictor representing its relative importance (RI; higher values represent superior relative importance) which refers to the total Akaike weight of all models containing that variable. When models included all reef habitat types, model selection was restricted to select either habitat or depth as collinearity precluded inclusion of both predictors. The best-fit model structure for each benthic response variable was retained, as well as the RI of each individual predictor.
Table 4.3. Explanatory models with their included predictors. *Depth was included in models (i) and (ii) focusing on outer reefs only, and habitat within models (i) and (ii) encompassing all reef habitat types. Best-fit model selection for all reef habitats was restricted to select either depth or habitat.
Model Comments Predictor set
(i) Spatial and physical available without reef
surveying depth*, habitat*, latitude/DHW, storm exposure
(ii) Local features available from reef
surveying browser biomass, depth*, excavator biomass, grazer biomass, habitat*, local human density, scraper biomass (iii) Best-fit from model selection
based on Akaike
Information Criterion (AIC)
selected from all possible predictors
To account for non-linear relationships between benthic assemblages and predictors, a generalised additive modelling (GAM) approach was adopted within the mgcv package in R (Wood, 2011). Spatial autocorrelation was checked using correlograms (ncf package in R – Bjornstad, 2016) for all response-model combinations, based on both response variables as well as model residuals. Spatial plots of the model residuals were also checked, with no correlations observed. All models (Table 4.3) were run separately for reefs exposed to different levels of local human impact, focused first on outer reefs (n = 62) and then for all reef habitat types (n = 182), where habitat replaced depth as a predictor in spatial and physical and local features models. To further explore differences in key predictors of benthic assemblages between reefs exposed to low and high local impacts, best-fit1 models tailored for low impact reefs were applied to high impact reefs (best-fit2).
To quantify the explanatory power of models, the adjusted-R2 (adj-R2) values were retained from each model. All data analyses were carried out in R version 3.1.1 (R Development Core Team, 2013), and graphs built using the ggplot2 package (Wickham, 2009).
73