Geomorphically controlled coral distribution in degraded shallow reefs of the Western Caribbean

Study region and geomorphic framework

The Mesoamerican Reef Tract (MART), located within the Western Caribbean marine ecoregion, stretches ∼1,000 km along the coasts of the Yucatan peninsula and northeast sector of Central America (Fig. 1; Spalding et al., 2007). It is composed of shallow detached fringing reefs that consist of two main geomorphic zones, a protected back-reef zone and an exposed reef-front zone separated by a crestline where waves break. These breakwater structures are limited to shallow water (<10 m deep) and are developed over and adjacent to a bedrock terrace which is veneered by coral grounds (Rodríguez-Martínez et al., 2011). Together the reef and terrace form a relatively consistent seascape over the inner fore-reef shelf (<15 m) throughout the MART (Blanchon et al., 2017). The shallow geomorphic framework of the MART is therefore consistent with the geomorphic zonation scheme first described from Grand Cayman by Blanchon & Jones (1995) and later summarized by Blanchon (2011), who both used slope breaks to identify reefal boundaries. We adopted this basic geomorphic scheme for all sectors of the MART including both Belize and Honduras.

In the northern Mexican sector of the MART, fringing reefs are relatively underdeveloped with limited reef-fronts that only extend to a depth of ∼6−8 m, thereby producing a more extensive bedrock terrace veneered by coral grounds (Blanchon et al., 2017). The edge of this terrace is marked by a slope break (mid-shelf break) and is commonly covered by a non-accretionary coral-ground community (Rodríguez-Martínez et al., 2011). In the Belize sector, fringing reefs transition into better developed Barriers and Atolls, and their reef-fronts extend into deeper water covering the bedrock terrace. At some sites, reef-front subzones like the spur-and-groove can extend all the way to the shelf edge in 25–30 m of water (James & Ginsburg, 1979; Burke, 1982). In the Honduras sector, reef development is restricted to the offshore Bay Island group: Roatan, Utila, Guanaja, and Cayos Cochinos. These reefs also have extensive spurs-and-groove zones extending to the shelf edge (Almada-Villela et al., 2003).

Following the geomorphic scheme of Blanchon (2011) we classify all sites based on retrospective satellite imagery from virtual globes such as Google Earth and ArcGIS online. A small number of sites (∼2%) have unclear images due to cloud cover and/or sea state, and geomorphic zoning is more difficult to delineate at the subzone level (for example, defining the boundary between the spur-and-grooves and non-tree bedrock terrace). In these few cases, we classify the sites as the most common type for that depth and zone based on field experience and expert opinion. Sites where reefs are less organized (∼17%) and lack a clear geomorphic structure, are simply classified as undifferentiated zones and assigned as either protected or exposed wave environments. Although uncommon, some of these ‘exposed’ zones can be found in semi-sheltered sites with complex coastal geography (e.g., the barrier reef fronts behind offshore atolls in the Belize sector).

Coral distribution in geomorphic zones

To document the spatial distribution of scleractinian coral within shallow geomorphic zones we use information from two databases, the Caribbean Reef Information System (CRIS, Barcolab.org) and the Healthy Reefs Initiative (www.healthyreefs.org), which contain qualitative data on coral cover estimates at site level. Data in these sources come from published literature, research projects, and monitoring programs. In all cases, coral cover information was obtained using two similar benthic survey methodologies (Line Intercept Transect- LIT and Point Intercept Transect-PIT) according to Atlantic and Gulf Rapid Reef Assessment Benthic Protocol (AGRAA, Lang et al., 2012 available at https://www.agrra.org/coral-reef-monitoring) and the Mesoamerican Barrier Reef System monitoring program (Almada-Villela et al., 2003). Each survey includes between six and ten transects (10 or 30 m in length) haphazardly deployed at the reef sites. From these databases we select only the most recent data from sites in a 1–10 m depth range that contain coral cover estimates at species level. After discarding sites outside these requirements (e.g., depth range, data quality at species level), we compile a subset of coral coverage from 95 survey sites between 2016 and 2018 (Table 1). This time period is chosen to avoid potential bias caused by the recent outbreak of the Stony Coral Tissue Loss Disease reported in the Mexican sector during 2018 (Alvarez-Filip et al., 2019).

Based on this subset, we calculate (i) the overall site live coral-cover data (LCC) and (ii) the relative coral cover of each species (i.e. the relative contribution of each species to the total coral cover). To avoid potential confusion stemming from using the term ‘coral cover’ to refer to the relative coral cover of each species, hereafter we refer to the latter as ‘coral species contribution’. In ecological studies, this is referred to as the sites taxonomic β-diversity and can be partitioned into two components: one reflecting species replacement between communities (turnover) and other reflecting the difference in the number of species among communities (nestedness-resultant; Baselga, 2012). However, coral-cover data are pooled and analyzed according to their geomorphic zonation, and not survey sites. Therefore, we compare overall coral-cover data of multiple sites as a comparative sampling unit, for example, pooling coral-cover from all back-reef sites (Fig. 2B). The qualitative β-diversity analysis entails partitioning the total dissimilarity in species composition (represented by the Sorensen index) into its turnover and nestedness components using 100 samples of 10 sites from each geomorphic zone. The analysis was conducted on Sorensen dissimilarity matrix using a betapart package available in R CRAN (Baselga et al., 2017).

Statistical analysis

The statistical approach we use for the comparative analysis of benthic parameters and their geomorphological settings is similar to the method described by Medina-Valmaseda et al. (2020) who analyzed the geomorphic patterns of coral species at a single site, Punta Maroma. The main goal in using multivariate analysis in that study, and in this one, is to minimize method-bias of benthic surveys from multivariate-based analyses of ecological, abundance-based variables, which include a narrow set of pre-treatments and transforming of the raw data (e.g., standardization). For the multivariate pre-treatment, we choose a medium level (square-root) of data transformation. Further transformation is followed by an ordination process where we construct the correspondent Bray–Curtis matrix of similarities. All subsequent multivariate statistical analyses and graphical outputs are performed and constructed using Plymouth Routines in Multivariate Ecological Research (Primer-e version 7.0.13, serial number 4901, Clarke & Gorley, 2015).

To compare overall coral cover and related coral species contribution to relative cover of the communities between geomorphic zones, we conduct an asymmetrical analyses of permutational variance (PERMANOVA) at both levels of aggregation for factor ‘site’ (nested within the geomorphic zones by wave exposure) and nested in wave exposure regimes on the basis of Bray–Curtis similarity measures of transformed square-root matrix of ecological data (Anderson, 2006). Including these permutations, procedures offer an alternative to Normal-based inferential methods as they are adaptable and require fewer assumptions, whereas the asymmetrical design deals with a different number of geomorphic zones by wave exposure, with three sheltered and four exposed. The experimental design consists of three factors (Fig. 2C): Factor A: wave exposure (fixed with a = 2 levels: sheltered and exposed), Factor B: Geomorphic zones (random, nested in wave exposure with five levels: irregular, coral-ground, spur & grooves, reef-front and back reef), and Factor C: Site (random, nested in geomorphic zones, wave exposure) with 95 levels. The test uses permutation of residuals under a reduced model and Type III (partial Square Sums) in 9,999 permutations.

To test the homogeneity of multivariate data dispersion in each case, we performed a non-parametric permutational analysis of multivariate dispersion (PERMDISP), along with pairwise comparisons of the Bray–Curtis matrix of similarities. PERMDISP is performed based on distances to centroids, with P-values (p(perm)) obtained from 9999 permutations (Anderson, Ellingsen & McArdle, 2006). To determine the species contributions to coral cover within each geomorphic zone we conduct a two-way similarity percentage analysis (SIMPER) for zones by wave exposure, based on Bray–Curtis similarity measures of transformed square-root matrix of abundance data, making a 70% cut-off for low contributions (Clarke & Warwick, 1994; Clarke et al., 2014). The results are presented through the metric MDS of bootstrapped averages for 95% region estimates. Bootstrap averages test iteratively resample each group by its geomorphic zone 200 times, creating a plot of 4 m dimensions multivariate effect. Such a bootstraps procedure allows better visualization of how the mean response of each data group varies with changes in predictor values if we were to collect another sample of observations from a different set of geomorphic zones (Fieberg et al., 2009). The arrangement of the number of sites by its geomorphic zone classification lacks sampling balance (Fig. 2B) and therefore, another advantage of this approach is its success in estimating parameters of a statistical distribution for a balanced number of copies from our data whereas preserving the original structure of its data set (Fieberg, Vitense & Johnson, 2020).

Extrinsic factors

For the purposes of this study, we consider extrinsic factors to be the diverse environmental variables that describe the abiotic context of coral communities including physicochemical variables and non-biotic disturbances, such as depth gradient, wave exposure and hurricane impacts. To evaluate the response of the biotic multivariate data to its corresponding multivariate environmental data we select information from both benthic surveys and cyclonic disturbances available for the MART. We use depth of benthic surveys and environmental information extracted from the open-source Physical Environments of the Caribbean Sea classification of Chollett et al. (2012) as geographic information system layers. For this, we extract data of physical disturbances including exposure to wind- and tropical-storm waves and the number of hurricane impacts using the QGIS (v. 3.14 Pi) Point Sampling Tool plugin, v.0.5.3 available at https://github.com/borysiasty/pointsamplingtool. Where the plugin failed because of no-value pixels, we manually selected the neighbouring major value. Furthermore, we check the distribution of all environmental variables for skewness and outliers and log-transform all data except the depth of each sampling site (following Clarke & Gorley, 2006). Finally, we normalize data to place each variable on the same dimensionless scale (for example, Clarke, Tweedley & Valesini, 2014). Environmental data are graphically represented through a Principal Coordinate Analysis (PCO) where vectors are the raw Pearson correlations of variables with the Principal Component Analysis (PCA) values (Torgerson, 1958; Gower, 1966).

To analyze and generate an exploratory hypothesis model on the spatial relationship between multivariate data clouds of the biotic response variables (coral species distribution) and that of regional environmental gradients, we test a distance-based linear model (DISTLM, Legendre & Anderson, 1999). The DISTLM test includes the Bray–Curtis matrices of previously square transformed data of coral species distribution and logarithmic Euclidean distance-based matrix of previously normalized environmental data in 9999 permutations within the factor Geomorphic zone for the environment. The test assumes non-linearity and additivity of the abiotic variables on the high-d community response. To visualize fitting models in multi-dimensional space we use the distance-based redundancy analysis (dbRDA, Legendre & Anderson, 1999; McArdle & Anderson, 2001) based on the Euclidean-distance matrix of environmental data.

Coral cover patterns in geomorphic zones

A total of 95 survey sites with 50 from Mexico, 25 from Belize and 20 sites from Honduras are included in this study. Sites encompass seven different types of geomorphic zones and subzones in both wave-exposed and protected environments along the regional latitudinal gradient and contain a total of 40 coral species. As shown in Fig. 3 results of β-diversity analysis show high levels of Total β-diversity (Fig. 3A, solid lines) for every geomorphic zone with high values and contribution of their turnover component (Fig. 3B, dashed lines) and low values for their nestedness-resultant component (Fig. 3B). Such results indicate that there is a low proportion of shared species between coral communities from adjacent geomorphic zones with a similar number of species. Although differences in taxonomic β-diversity values are clear between coral communities from sheltered back reef and exposed coral-ground zones, these are diffuse and overlapping for the communities from exposed reef front and spur & groove zones. The lower value of the nestedness-resultant fraction of Sorensen dissimilarity belongs to coral communities from sheltered back reef zones, whereas those from the irregular zones of wave-sheltered environments and wave-exposed spur & groove and coral-ground zones show almost identical patterns, indicating a confluence in these components for the exposed zones (Fig. 3B). Under the setting conditions of the β-diversity test, more sites were requested for sampling in the available datasets from two of the seven geomorphic zones. Results of the test are therefore unavailable for wave-exposed irregular zones (6 sampling sites), and wave-sheltered reef fronts (2 sampling sites).

The results of the PERMANOVA test indicates that there are no significant differences in overall coral cover between wave sheltered and wave-exposed geomorphic zones (P(perm) = 0.39; Fig. 4A; Data S1), although when species identity is included in the relative composition analysis, there are significant differences between geomorphic zones (Fig. 4B). The results of the asymmetrical PERMANOVA test including the three terms ‘wave exposure’, ‘Geo zone (wave exposure)’ and ‘site (Geo zone (wave exposure))’ show that there are marginal differences among the wave exposure environments in assemblage structure (P(perm) = 0.06; P(MC) = 0.07) whereas there is a significant variability among assemblage structure from the geomorphic zones (P(perm) and P(MC) < 0.01, Fig. 4A; Data S2). The estimates of components of variation (S) for the three orthogonal factors are wave exposure: 115.82; Geo zone nested in wave exposure: 214.98; and site nested in Geo zone nested in wave exposure: 875.33. The site factor therefore contributes most to differences in species contribution, followed by the geomorphic zone, with wave exposure having the least control.

A posteriori pairwise PERMDISP test at the LCC level shows homogeneity of variance in all groups (F: 1.49, df1: 6, df2: 84; p(perm): 0.3422; Data S3), whereas the species contribution level shows heterogeneity in data variance involving pairwise tests for single or combined irregular and exposed zones (F: 13.23, df1: 6, df2: 543; p(perm): <0.01; Data S4). Other cases, including the irregular and combined sheltered/exposed geomorphic zones, show homogeneity in data variances (F: 13.23, df1: 6, df2: 543; p(perm): <0.01; Fig. 4B). In total 61% of possible pairwise tests show homogeneity in data variance and do not show a particular trend under the current PERMANOVA design. This indicates that there are differences in coral-species distribution between geomorphic zones (PERMANOVA), and in the variance of data distribution (PERMDISP) in several zones within the same environments (Fig. 4, Data S4).

The contribution of individual species to observed differences is assessed using the SIMPER test (Data S5). Results show that in exposed geomorphic zones three species overlap: Agaricia agaricites, Porites astreoides and Siderastrea siderea (Average similarity: 30.0) with relative contributions of 32.2%, 24.3% and 15.1% respectively. In sheltered zones another trio of species overlap, P. astreoides, Orbicella annularis and Agaricia agaricites (Average similarity 25.7) with relative contributions of 49.0%, 18.5% and 8.4% respectively (Data S1 to Data S5). Interestingly, an almost identical group of species accounts for the 81.2% average dissimilarity between samples. In terms of dissimilarity, only four species (S. siderea, Porites astreoides, O. annularis and A. agaricites), which contribute ∼10%, form roughly 53% of the differences between samples. Overall, across all wave exposure zones, the largest average similarity of 34.9% corresponds to the exposed Coral-ground zone, and the lowest 18.2% to the sheltered Back-reef zone.

Coral species patterns and extrinsic factors

Average species richness (S, the number of species identified to lowest possible taxonomic level) increase with depth being the lowest in sheltered back-reef (2.1) and the highest in deeper coral-ground zone (4.9). Also species diversity as measured by the Shannon–Weiner index (H′) follows the same pattern, increasing seaward from the back- reef (0.9) towards coral-ground zone (1.4; Table 1).

On a regional scale, all sites are exposed to a similar wave-exposure regime regardless of geomorphic zone, whereas hurricane impact varies between sites (Fig. 5). When analyzed by geomorphic zone and wave exposure, the average site depth varies from 1.8 ± 2.9 m in sheltered reef-front zones, to 11.2 ± 4.2 m in exposed coral-ground zones. Chronic stress derived from wave-exposure for these sites varies between 5.9 ± 0.9 and 7.4 ± 0.2 J m⁻³. Whereas acute mechanical stress derived from hurricane frequency in the last 157 years (1851–2008) varies between 11 ± 2 years and 20 ± 2.8 years (Table 1). Regional patterns of three selected environmental data are presented in the Principal Coordinates Analysis (PCO) of Fig. 5. It shows that the first two axes explain 81.6% of the internal variability of these three environmental variables, indicating that the 2-d ordination is likely to have captured the majority of the patterns of multidimensional data variation. Accordingly, the depth of the reef site and hurricane impact explain 42.2% of variations of the PCO1 axis, and wave energy exposure alone explains 39.5% of the PCO2 axis. In addition, depth and hurricane impact appear to act inversely over coral species patterns (Fig. 5A).

Results from the DISTLM marginal tests (Data S6) show how much can be explained by each abiotic variable alone, ignoring other variables. Tests show that the depth of reef sites exerts the highest influence explaining 5.5% of the overall variability in coral species contribution (p < 0.01, Data S4A, Fig. 5B). Wave exposure and the number of major hurricane impacts when each is considered alone explain roughly 1% to 1.5% (p < 0.01, Data S4A, Fig. 5B) of variability in species contribution to each geomorphic zone. These results indicate that spatial patterns in coral species contribution by geomorphic zones are still responsive to these three environmental variables.