Multi-scale processes drive benthic community structure in upwelling-affected coral reefs

Environmental processes acting at multiple spatial scales control benthic community structures in coral reefs. However, the contribution of local factors (e.g., substrate availability and water clarity) vs. non-local oceanographic processes (e.g. upwelling events) in these highly complex systems is poorly understood. We therefore investigated the relative contribution of local and non-local environmental factors on the structure of benthic groups and specifically on coral assemblages in the upwelling-affected Tayrona National Natural Park (TNNP, Colombian Caribbean). Coral-dominated communities were monitored along with key environmental parameters at water current-exposed and -sheltered sites in four consecutive bays. Regression tree analyses revealed that environmental parameters explained 59.1% of the variation within the major benthic groups and 36.1% within coral assemblages. Findings also showed recurring patterns in community structures at sites with similar exposure across bays. We suggest that benthic community composition in TNNP is primarily driven by 1) wave exposure, followed by 2) temporal changes in nutrient availability governing the structure of benthic groups, and 3) local bay-specific differences controlling the zonation of benthic groups and coral assemblages. This study highlights the existence of complex hierarchical levels of local and non-local environmental factors acting on reef communities and stresses the importance of considering processes operating at multiple spatial scales in future studies on coral reef community structure and resilience.


INTRODUCTION
Large-scale abiotic parameters including global water temperatures and latitude-dependent light availability control the distribution of tropical coral reefs (Veron, 1995;Kleypas et al., 1999). They are therefore generally restricted to shallow, warm and clear waters at latitudes between 30 • N and 30 • S (Kleypas et al., 1999;Veron and Stafford-Smith, 2000). In contrast, regional processes such as ocean circulations and upwelling events that operate on smaller scales account for coral reef distribution within their global limits and for the heterogeneity of reef community structures by controlling reef species composition (Hubbard, 1996;Karlson and Cornell, 1999;Pandolfi, 2002). In addition, local factors such as substrate availability and sedimentation rates can vary at lower spatial scales and thus determine the variability in organism abundances, thereby contributing to the small-scale heterogeneity in community structures (Geister, 1977;Pandolfi, 2002). However, some factors interact with each other across spatial scales and are thus difficult to classify. Wave exposure, for example, may change at scales of a few meters, depending on local geomorphologies such as coastline and seafloor structures, but is largely controlled by prevailing wave and current regimes (Hamner and Wolanski, 1988;Ekebom et al., 2003;Monismith, 2007). Therefore, wave exposure essentially is a regionally-driven local factor. As a consequence, the precise estimation of the contribution of local vs. regional factors to reef community structure remains a challenge (Done et al., 1991;Karlson and Cornell, 1999;Edmunds, 2013). This situation becomes even more complicated due to the fact that abiotic factors may strongly affect the supply-side ecology (e.g., fecundity and dispersal ability) of sessile organisms and therefore their community structure (Lewin, 1986;Underwood and Fairweather, 1989;Hughes et al., 1999;Jones et al., 2009).
Moreover, a proper statistical assessment of the effects of local and non-local factors on coral reef community structures requires systems that ideally comprise a number of comparable units that are affected by regional factors but also show specific local differences. Unfortunately, only few systems that fulfill these criteria are available for studying reef benthic community structure: one example are the reef communities of the Similan Islands (Thailand), which are strongly affected by pronounced short-term variations in water temperature, pH, and nutrient availability due to large-amplitude internal waves (LAIW, Schmidt et al., 2012). Although the reef communities of the islands follow similar ecological functions, the environmental variations induce spatial differences of benthic community structures between LAIW-exposed and -sheltered sites (Schmidt et al., 2012). Another potential system are the Pearl Islands at the Pacific coast of Panamá, where upwelling-induced gradients in water temperatures account for similar patterns in benthic reef communities with limited coral growth at the upwelling exposed sides of the islands (Glynn and Stewart, 1973).
Here, we propose the Tayrona National Natural Park (TNNP) at the Colombian Caribbean coast as a promising system for the investigation of the effects of local and regional environmental factors on coral-dominated benthic communities. Within TNNP, several ecologically similar bays are located in close proximity to each other (Figure 1) that harbor diverse coastal ecosystems including highly structured coral-dominated benthic communities and coral reefs (Garzón-Ferreira and Cano, 1991;Díaz and Acero, 2003;Garzón-Ferreira and Díaz, 2003). During the past two decades, benthic communities remained relatively stable with no major disturbance events affecting the area (Rodríguez-Ramírez et al., 2010): coral bleaching was negligible (Navas-Camacho et al., 2010) and only one hurricane led to minor reductions in coral cover in 1999 (Rodríguez-Ramírez et al., 2010). However, despite this relative lack of disturbances, distinct intra-and inter-bay differences in benthic community structure have been reported in TNNP (Garzón-Ferreira and Vega-Sequeda et al., 2008).
One of the major processes that affect the bays is seasonal upwelling (Antonius, 1972;Diaz-Pulido and Garzón-Ferreira, 2002;Eidens et al., 2014), leading to significant changes in environmental parameters such as water temperature, salinity, and nutrient availability (Bayraktarov et al., 2014a). Upwelling-related environmental changes in these parameters may, on the one hand, adversely affect coral growth (Antonius, 1972) and induce macroalgal blooms (Diaz-Pulido and Garzón-Ferreira, 2002). On the other hand, they may also cause a decrease in bleaching susceptibility of corals (Chollett et al., 2010) and facilitate recovery from coral bleaching (Bayraktarov et al., 2013). In addition, environmental changes also influence primary production of major benthic photoautotrophs (Eidens et al., 2014), thereby affecting benthic community composition. Other factors that potentially influence benthic communities in the area include wave exposure, substrate availability, and changes in light intensity (Antonius, 1972;Erhardt, 1975;Werding and Sánchez, 1989;Garzón-Ferreira, 1998;Vega-Sequeda et al., 2008;Bayraktarov et al., 2014b).
Despite several studies on potential effects of individual environmental parameters on benthic communities in TNNP, it is still not fully understood which factors contribute to the present patterns in benthic community composition and at which scales these factors operate.
The major goal of this paper therefore was to use multivariate models in order to identify the drivers of benthic community structures across spatial scales in TNNP coral reefs. The study consisted of the following working steps.
(1) Two contrasting sites (water current-exposed and -sheltered) in four adjacent bays (Figure 1) were selected, based on preliminary investigations (conducted in 2010) and on former benthic community surveys (Garzón-Ferreira and Cano, 1991). Wave exposure was quantified for each study site.
(2) The community structures of major benthic groups/substrates including scleractinian corals, algal turfs, frondose macroalgae, crustose coralline algae, sponges, and sand (hereafter referred to as benthic groups) were monitored together with important environmental parameters (i.e., temperature, salinity, nutrient availability, and water clarity) during upwelling and non-upwelling periods.
(3) Scleractinian corals were identified to the species level (hereafter called coral assemblages). (4) Generalized linear mixed-models (GLMMs) were used to classify the spatial scales at which the environmental parameters operate. (5) In a final step, multivariate regression tree (MRT) and indicator species analyses (ISA) were used to identify (a) the processes and parameters driving benthic community structure together with their individual contribution, and (b) the benthic groups being most affected by these drivers.
Our working hypothesis was that wave exposure causes, on the one hand, significant differences in community structures within bays (i.e., exposed vs. sheltered sites) but, on the other hand, similarities within groups of exposed vs. sheltered sites across bays.

STUDY AREA
The TNNP is located at the Caribbean coast of Colombia near the city of Santa Marta. Due to the proximity to the Sierra Nevada de Santa Marta, the world's highest coastal mountain range, the coastline is heterogeneous and interspersed with numerous bays and inlets (Figure 1). The climate in the region is strongly influenced by the seasonal Caribbean Low-Level Jet of northeast trade winds (Andrade and Barton, 2005). Wind incidence is highest during the dry season from December to April (<20% of the annual rainfall) and generates a coastal upwelling through Ekman transport of surface waters away from the coast (Andrade and Barton, 2005). During this period (hereafter referred to as upwelling season), strong wind-induced water currents and wave action prevail, mean monthly water temperatures decrease to around 25 • C with minimal water temperatures being as low as 20 • C, and salinity as well as nutrient availability increase (Bayraktarov et al., 2014a). During the major rainy season from September to November (hereafter referred to as non-upwelling season), trade winds cease and surface water temperatures increase (monthly means around 28 • C; maximum values greater than 30 • C, Salzwedel and Müller, 1983;Bayraktarov et al., 2014b). During this period, increased precipitation and terrestrial run-off may result in higher turbidity and lower salinity (Salzwedel and Müller, 1983;Diaz-Pulido and Garzón-Ferreira, 2002).
The study was conducted in four adjacent bays of TNNP (Chengue, Gayraca, Neguanje, and Cinto; Figure 1). These bays are comparable regarding their environmental settings (e.g.,

FIGURE 1 | Location of investigated bays in the Tayrona National Natural
Park at the Caribbean coast of Colombia. White circles indicate study sites. Wind rose shows wind speed and direction for the area (QuickSCAT data, 1999(QuickSCAT data, -2008. Abbreviations of study sites: Ch S , Chengue sheltered; Ch E , Chengue exposed; Ci S , Cinto sheltered; Ci E , Cinto exposed; Ga S , Gayraca sheltered; Ga E , Gayraca exposed; Ne S , Neguanje sheltered; Ne E , Neguanje exposed. Data source: digital elevation model from ASTER GDEM by MEDI (Japan) and NASA (USA), bathymetry from Centro de Investigaciones Oceanográficas e Hidrográficas del Caribe (Colombia), coral formations from Garzón-Ferreira and Cano (1991) and INVEMAR (2012).  Chollett and Mumby (2012) taking into account fetch length and wind speed and direction. Anthropogenic impact was quantified using selected categories of the IUCN threat classification scheme (IUCN, 2001) and the scoring system provided by GIWA (2001) with scores ranging from 0 (no known impact) to 3 (severe impact).
Prevailing waves and surface water currents, predominantly coming from NE, are strongly affected by local topographies.
Thus, similar orientations of the bays result in comparable patterns of wave exposure (Werding and Sánchez, 1989). Whereas the outer western parts of the bays are strongly exposed, the eastern parts are generally sheltered.
One study location with comparable coral communities at the exposed western and the sheltered eastern side of each of the four bays was chosen, resulting in a total of eight study locations (Figure 1). For each study location, wave exposure (in J·m −3 ) was calculated using the method of Chollett and Mumby (2012), taking into account the distance of water over which the wind blows together with wind direction and speed. Briefly, fetch (i.e., the openness) was quantified by tracing a line in 48 compass directions from each of the study locations across the sea until land (as defined by the Global Administrative Areas database) was encountered. Daily wind speed and direction over a 10-year period (i.e., 1999-2008) and a spatial resolution of ∼12.5 km were taken from the QuickSCAT Project (http://podaac.jpl.nasa. gov/dataset/QSCAT_LEVEL_2B_OWV_COMP_12). Fetch and wave exposure were calculated in R 3.0.3 (R Core Team, 2014) including the packages sp 1.0-15 (Pebesma and Bivand, 2005) and ncdf 1.6.6 (Pierce, 2011).
As community structures in TNNP may vary over time depending on season (i.e., upwelling vs. non-upwelling, Diaz-Pulido and Garzón-Ferreira, 2002), benthic community assessments took place during the non-upwelling season in November 2011 and at the end of two upwelling seasons in April 2011 and 2012 to account for natural fluctuations in environmental parameters and/or benthic community structure in the study area. Note that we conducted two post-upwelling surveys as the effect of upwelling events on benthic communities may vary from season to season (Diaz-Pulido and Garzón-Ferreira, 2002;Eidens et al., 2014).

Environmental parameter assessment
Water temperatures at all eight study sites were measured continuously in situ every 5 min by calibrated HOBO TidBit v2 temperature loggers or Pendant temperature loggers (Onset Computer Corp., Bourne, USA) attached to the reef structure at a water depth of 10 m. For subsequent analyses, we used mean monthly water temperatures during the time when benthic monitoring took place. Salinity was measured at each study site using a portable meter (HQ40d, Hach, Loveland, USA) equipped with a 4-pole conductivity probe (CDC401, Hach, Loveland, USA). At all sites, water samples (n = 3) for inorganic nutrient determination [nitrate (NO − 3 ), nitrite (NO − 2 ), and soluble reactive phosphorus (SRP) mainly in the form of orthophosphate] were collected from the water column at ca. 1 m above the seafloor during each monitoring campaign. Samples were cooled on ice during transportation, immediately filtered in the laboratory (glass micro fiber filters, 0.7 µm particle retention, VWR International, Radnor, USA) and frozen until further analyses. Nutrient concentrations were quantified spectrophotometrically according to Garay et al. (2003). Water clarity was evaluated monthly by Secchi discs at each site (n = 4), and pooled means for upwelling/non-upwelling seasons were used. Data on environmental parameters and high-resolution temperature is available at pangaea.de (Bayraktarov et al., 2014c,d). We visualized the multivariate variability of z-standardized water parameters by a principal component analysis (PCA) using the R package vegan 2.0.10 (Oksanen et al., 2012). Anthropogenic impact within investigated bays was evaluated using selected categories of the IUCN threat classification scheme (IUCN, 2001; Table 1) and the scoring system provided by GIWA (2001) with scores ranging from 0 (no known impact) to 3 (severe impact).

Benthic community assessment
In order to monitor benthic communities, we established three fixed 50 m transects at the 10 m isobath in each of the eight study sites. Transects were separated by gaps >5 m to ensure independence among samples. Start and end points of each transect were marked with underwater buoys, and transects labeled with cable ties on the seafloor in 3-5 m intervals. Benthic community structures along the established transects were determined using the line point intercept method (modified from Hodgson et al., 2004). Benthic groups were monitored at 0.5 m intervals directly below the measurement point and were categorized into hermatypic corals (including reef building hydrozoans of the genus Millepora), frondose macroalgae, algal turfs (sensu Steneck, 1988), crustose coralline algae (CCA), sponges, other biota, sand, and abiotic substrate. For the detailed assessment of coral assemblages, the same methodology was used as for benthic groups, but taxa were identified at the species level.
To display a low-dimensional figure of community composition of benthic groups and coral assemblages, non-metric multidimensional scaling (NMDS, Kruskal and Wish, 1978) was utilized. We used the Morisita-Horn dissimilarity index on log(x + 1) transformed benthic cover to perform a three-dimensional NMDS with the R package vegan 2.0.10.

ASSESSING SPATIAL SCALES OF QUANTIFIED ENVIRONMENTAL PARAMETERS
GLMMs (Bolker et al., 2009) with fixed and random effects were used to assign spatial scales to quantified water parameters. We used the predictor variables study period, bay, and exposure within bays (i.e., wave-exposed and -sheltered side of each of the four bays) as fixed effects to explain the variance in both a multivariate GLMM with all water parameters (z-standardized), and several univariate GLMMs with water temperature, salinity, nutrient availability, and clarity as individual response variable. Because of repeated measurements at the same sites in the same four bays, we nested sites within bays. This random effect of the model structure accounts for the statistical non-independence of our observations (Bolker et al., 2009). Since appropriate GLMM test statistics are under debate (Bolker et al., 2009), here we applied the Markov chain Monte Carlo technique as implemented in the R package MCMCglmm 2.19 (Hadfield, 2010). In order to infer the significance of model predictors from the posterior samples, we used the default weakly-informative prior, a burn-in of 10,000, and a thin of 25 applied to 60,000 Markov chain iterations. Explained variance of the three predictors was calculated according to Nakagawa and Schielzeth (2013).
We applied the following operational criteria to estimate a posteriori whether a parameter is driven by regional or local processes: a high percentage of explained variance among sites and bays classifies a factor as local, strong recurring intra-bay variation indicate local factors driven by larger-scale processes that affect all bays similarly. If spatial differences between investigated sites and bays explain none or only a minor proportion of variance in the data, we classified the variable as non-local factor.

IDENTIFYING ENVIRONMENTAL PARAMETERS STRUCTURING BENTHIC COMMUNITIES
MRT analyses (De'ath, 2002) were used to identify driving environmental factors and their individual contribution to benthic community structure. This method performs hierarchical dichotomous clusterings of community data by selecting environmental parameters (including water temperature, salinity, NO − 3 , NO − 2 , SRP, wave exposure, water clarity, and other nonquantified bay-specific factors) that maximize the homogeneity within groups of line transects. Accordingly, these clusters are characterized by both a homogenous assemblage structure and similar ecological parameters. MRTs are not based on traditional significance testing but on 10-fold cross-validation (CV) in order to determine the number of dichotomous splits and the importance of predictor variables (De'ath, 2002). We used the R package mvpart 1.6-0 (Therneau et al., 2012) to perform the analyses with 1,000 cross-validation runs on Hellinger transformed benthic group and coral assemblage data. Since co-linearity among variables may impair the selection of regional and local factors during MRT analyses, Pearson's correlation analysis was performed in R to test for correlations among water parameters. In case of correlation coefficients >0.8, the ecological parameter offering an already known effect on benthic communities was retained in the analyses. This standard procedure in ecological modeling reduces the amount of explaining factors and ensures a balanced design without overrepresentation of correlated factors (Zuur et al., 2010). Additionally, ISA (Dufrene and Legendre, 1997) was performed using the R package MVPARTwrap 0.1-9.2 (Ouellette and Legendre, 2013) in order to identify benthic groups and coral assemblages associated with splits derived from the MRT analyses. This method aims at detecting species that represent distinct ecological settings and indicate site-specific community types.

ENVIRONMENTAL PARAMETER ASSESSMENT
The first two axes of the PCA, which was conducted to visualize the multivariate variability of z-standardized water parameters, explained 75.1% of the variation (Figure 2; for exact values at each study location see Supplementary Material Tables S1-S4). The first axis could mainly be related to changes between study periods with strong variation in salinity, NO − 3 and NO − 2 availability, and water temperature. Whereas salinity as well as NO − 3 and NO − 2 availability were higher during upwelling events, water temperatures were higher during non-upwelling.
The GLMM analyses, performed to assign a posteriori spatial scales to water parameters, revealed differences in means of all investigated water parameters but no significant intra-bay differences for any factor ( Table 2). In contrast, inter-bay differences could explain variation in all water parameters except for SRP and water temperature. Differences between bays were the main source of variation in water clarity (54.5%) and explained variation in NO − 3 and NO − 2 availability to some extent (9.7 and 8.9%, respectively), but only 0.8% of variation in salinity. Fluctuations between study periods were the main source of variation in www.frontiersin.org February 2015 | Volume 2 | Article 2 | 5

FIGURE 3 | Graphical representation of (A) non-metric multidimensional scaling (NMDS) and (B) multivariate regression tree analysis (MRT) for major benthic groups in Tayrona National Natural Park. (A) Variation in
benthic community structure shown for the first two axes from a NMDS and the relative contribution of each benthic group (arrows). Colors indicate groupings according to the coloring of MRT leafs. (B) MRT and associated physical and environmental variable thresholds related to the separation of major benthic groups separating individual transects according to wave exposure (in J·m −3 ), nitrite and soluble reactive phosphate availabilities (in µmol·L −1 ), and bay-specific factors. R 2 is the variance among the communities explained by the respective environmental parameter. NO 2 , nitrite; SRP, soluble reactive phosphate; Ch S , Chengue sheltered; Ch E , Chengue exposed; Ci S , Cinto sheltered; Ci E , Cinto exposed; Ga S , Gayraca sheltered; Ga E , Gayraca exposed; Ne S , Neguanje sheltered; Ne E , Neguanje exposed.
all investigated parameters except for water clarity (water temperature 99.8%, salinity 83.6%, NO − 3 availability 63.8%, NO − 2 availability 53.4%, SRP availability 13.7%). We therefore classified water clarity as a mainly local factor and all other water parameters as mainly non-local factors.

BENTHIC COMMUNITY ASSESSMENT
Our transect data revealed that scleractinian corals, turf algae, CCA, frondose macroalgae, and sponges were the dominant benthic groups at the monitored sites and, together with sand, amounted to 96.2 ± 0.6% of total benthic coverage (see Tables  S6-S8). Among these groups, corals and algal turfs accounted for the highest coverage (24.8 ± 1.2 and 23.4 ± 1.1% mean benthic coverage, respectively). In our community composition analysis, the first NMDS axis (NMDS 1, Figure 3A) differentiated transects along a gradient from high CCA (up to 45% per transect), coral (up to 52%), and sponge coverage (up to 16%), to high sand coverage (up to 67%). Within the second NMDS axis, transects were separated along a gradient from high macroalgae coverage (up to 50%) to high coverage of turf algae (up to 48%; NMDS 2).
As for the coral assemblages, a total of 25 species was identified during the study with a mean species richness of 12.3 species per study site. Minimum and maximum numbers of coral species were found in Chengue Bay with 10 and 17 species at the exposed and sheltered site, respectively (Table S9). The first NMDS axis (NMDS 1; Figure 4A) separated transects with high abundances of Montastraea cavernosa (up to 16% benthic cover) and Siderastrea siderea (up to 7%) from transects with high coverage of the Orbicella annularis complex (sensu Budd et al., 2012;up to 15%) and Colpophyllia natans (up to 17%). The second NMDS axis differentiated transects dominated by Pseudodiploria strigosa (sensu Budd et al., 2012;up to 30%) from transects with higher abundances of Millepora complanata (NMDS 2).

EFFECTS OF REGIONAL AND LOCAL ENVIRONMENTAL PARAMETERS ON BENTHIC COMMUNITY STRUCTURES
Prior to the MRT analyses, a Pearson's correlation analysis was conducted to test for correlations among quantified water parameters. It revealed a strong positive correlation between NO − 3 and NO − 2 concentrations (r Pearson = 0.8, Table S10) and a strong negative correlation between NO − 3 availability and water temperature (r Pearson = −0.81). We therefore excluded NO − 3 availability from the subsequent multivariate analyses.
For major benthic groups, a MRT analysis (CV-error = 0.65) revealed that wave exposure within bays, NO − 2 availability, SRP availability, and bay-specific factors together explained 59% of variation in benthic group structures ( Figure 3B). As indicated by the hierarchical order and branch length of the MRT, wave exposure was the main factor explaining 21.5% of the variation alone (split 1; Figure 3B), separating all locations at the exposed sides from the ones at the sheltered sides of the bays. The ISA revealed that this dichotomous split was mainly associated with significant differences in coverage of sand, sponges, and CCA (indicator values 0.67, 0.66, and 0.6, respectively, p ≤ 0.001 for all values). Whereas more sand coverage was present at sheltered sites, CCA and sponge coverage were higher at the exposed sites. Benthic groups at sheltered and exposed sites were further separated by differences in NO − 2 availability (splits 2 and 3; Figure 3B) that together accounted for 17.3% of variance in composition. These splits were, among others, associated with relatively higher macroalgal coverage at high NO − 2 concentrations (indicator value 0.7 and 0.79 for splits 2 and 3, respectively; p ≤ 0.001 for both values). In general, nutrient concentrations (NO − 2 and SRP availabilities) together explained 24.8% of benthic group structure (splits 2, 3, 5, and 15; Figure 3B). Local bayspecific factors accounted for 12.7% of variation (splits 4, 6, 7, 11, and 14; Figure 3A), and four out of these splits (4, 6, 7, and 11) were associated with differences in general coral abundance (indicator values 0.6, 0.56, 0.58, and 0.61, respectively, p ≤ 0.01 for all values).
Regarding coral assemblages, wave exposure and bay-specific factors together explained 36.1% of the variation in coral species composition (CV-error 0.64; Figure 4B). As for benthic groups, wave exposure was the main parameter to separate coral assemblages and accounted for 24.9% of variance (splits 1 and 4; Figure 4B). ISA revealed that the first split was mainly associated with patterns in abundance of P. strigosa, S. siderea, and the O. annularis complex (indicator values 0.6, 0.58, and 0.76, respectively, p = 0.001 for all values). Whereas P. strigosa and S. siderea generally exhibited higher abundances at more exposed sites, the O. annularis complex was more abundant at more sheltered sites (Figure 4A). Local bay-specific factors further discriminated coral assemblages from sheltered vs. exposed sites and together explained 11.1% of variation (splits 2 and 3; Figure 4B). At the sheltered sites, coral assemblages in Neguanje were separated from the sheltered sites in Chengue and Gayraca ( Figure 4A; split 2, Figure 4B), associated with higher coverage of C. natans and the O. annularis complex in Chengue and Gayraca (indicator value 0.77 and 0.63, respectively, p = 0.01 for both values). At the more exposed sites, the split between Neguanje exposed ( Figure 4A) and the other bays (split 3, Figure 4B) was mainly associated with higher abundances of S. siderea, Acropora palmata, C. natans, and Agaricia humilis in Neguanje (indicator values 0.7, 0.67, 0.63, and 0.59, respectively, p ≤ 0.01 for all values). Finally, wave exposure separated coral assemblages from the exposed sites in Chengue, Gayraca, and Cinto from Cinto sheltered ( Figure 4A; split 6, Figure 4B). This split was mainly associated with higher coverage of C. natans and the O. annularis complex at Cinto sheltered (indicator values 0.77 for both, p ≤ 0.001 for both values).

DISCUSSION
Our study revealed that environmental parameters operating on both local (pronounced differences between sites/bays in TNNP) and non-local scales (no/marginal differences between investigated sites/bays) could explain a considerable portion of variation in community structure of major benthic groups, whereas variation in coral assemblages could be mainly linked to local factors. The MRT analyses suggested that distinct differences in major benthic groups and coral assemblages between exposed and sheltered sites are present in all investigated bays and that community structure is generally similar at exposed as well as at sheltered sites (Figures 3B, 4B). Fluctuations in nutrient availabilities, largely occurring at non-local scales, further accounted for patterns in major benthic groups at all study sites, whereas local bay-specific factors explained a similar amount of variation in both, benthic groups and coral assemblages. Our results largely support the working hypothesis that wave exposure is the main driver of benthic community structure in TNNP. In the following, we and bay-specific factors. R 2 is the variance among the communities explained by the respective environmental parameter. Ch S , Chengue sheltered; Ch E , Chengue exposed; Ci S , Cinto sheltered; Ci E , Cinto exposed; Ga S , Gayraca sheltered; Ga E , Gayraca exposed; Ne S , Neguanje sheltered; Ne E , Neguanje exposed.

www.frontiersin.org
February 2015 | Volume 2 | Article 2 | 7 discuss these findings in more detail with respect to the specific effects of the investigated environmental parameters on benthic community structures.

EFFECTS OF LOCAL AND NON-LOCAL ENVIRONMENTAL PARAMETERS ON BENTHIC COMMUNITY STRUCTURES
Strong fluctuations in water temperatures, salinity, and nutrient availabilities were present at all study sites ( Table 2), being mainly associated with coastal upwelling events (Garzón-Ferreira, 1998;Diaz-Pulido and Garzón-Ferreira, 2002;Bayraktarov et al., 2014a). However, none of the quantified water parameters showed significant intra-bay differences ( Table 2; see also Bayraktarov et al., 2014b). On the contrary, similar topographic conditions among the bays in TNNP (Table 1, Figure 1) cause comparable gradients in wave exposure within the bays ( Table 1; Werding and Sánchez, 1989;Garzón-Ferreira and Cano, 1991). These gradients correlate with differences in benthic community structure between exposed and sheltered as suggested in our working hypothesis and presumed for coral assemblages in TNNP by Werding and Sánchez (1989). Wave exposure and water currents were previously shown to control the zonation of characteristic coral reef structures (Geister, 1977;Done, 1983) by influencing many ecological aspects of coral reefs such as water quality, sedimentation patterns, nutrient uptake, primary production, recruitment, larval dispersal, and bleaching patterns (Hamner and Wolanski, 1988;Andrews and Pickard, 1990;Tribble et al., 1994;Nakamura and Van Woesik, 2001;Nakamura et al., 2003;Monismith, 2007;Mass et al., 2010). ISA analyses revealed that sand, CCA, and sponges were the main benthic groups affected by gradients of wave exposure within the bays with higher sand abundance at the sheltered sites and higher CCA and sponge coverage at the exposed sites. Sediment transport and accumulation are generally correlated with hydrodynamics, as shown for a Hawaiian fringing reef where wave energy and water currents control sediment resuspension and transport (Ogston et al., 2004). Differences in sand coverage between sheltered and exposed sites in TNNP bays are therefore likely due to generally higher wave energy and water current velocities at exposed sites (Bayraktarov et al., 2014b), resuspending sediments and transporting resuspended material to sheltered or deeper parts where it accumulates . Our findings are in accordance with former studies, stating that sand abundance is inversely correlated with wave exposure (Bradbury and Young, 1981;Graus and Macintyre, 1989;Roberts et al., 1992). On the contrary, sponges and CCA are relatively resistant against hydrodynamic stress (Palumbi, 1986;Fabricius and De'ath, 2001) and are thus generally more abundant at reef sites exposed to high hydrodynamic forcing (Littler and Littler, 1984;Steneck, 1986).
Similar to the benthic groups, coral assemblages also differed significantly according to the degree of wave exposure (Figure 4B), which is in accordance with former studies, stating that coral community composition is commonly dictated by wave exposure (Geister, 1977;Done, 1983). Main reef-building corals that occur in TNNP (i.e., the O. annularis complex and C. natans) are generally restricted to the sheltered parts of the bays. Although other corals are abundant in exposed areas (i.e., P. strigosa and M. cavernosa), it seems that the main reef-building corals are not able to tolerate these turbulent environments due to recruitment failure, as suggested for the O. annularis complex by Chollett and Mumby (2012). The lack of main reef-building corals is an important factor that influences framework-building potential and may therefore largely account for the scarcity of reef frameworks at the exposed sites, together with temperature-induced decreases in coral growth (Glynn and Stewart, 1973) and possible seasonal abrasion during upwelling (Geister, 1977). Our assumption is supported by Grigg (1983) and Benzoni et al. (2003), suggesting that the development of reef frameworks in the Indo-Pacific is generally restricted to sites that are sheltered from high wave energy and/or upwelling.
As wave exposure is driven by current-and wave regimes that interact with local geomorphologies, we suggest that regional wave patterns account for the recurring patterns in community structure with similar differences in community structure present in the investigated bays and bay-independent similarities among exposed and sheltered sites. Comparable patterns in benthic community structure at the Pearl Islands (Panamá) and Similan Islands (Thailand) could also be explained by the exposure to large scale processes (coastal upwelling and LAIW, respectively; Glynn and Stewart, 1973;Schmidt et al., 2012). This highlights the effect of non-local processes in structuring benthic communities on small spatial scales.
The second level of separations in the MRT analysis of benthic groups (splits 2 and 3; Figure 3B) as well as subsequent splits were associated with upwelling-related fluctuations in NO − 2 availability and co-correlated parameters (water temperature, salinity, and NO − 3 availability; Table S10). These mainly non-local parameters thus constitute the second-strongest driver of benthic group composition in TNNP, primarily causing changes in algal assemblages with macroalgal blooms during upwelling events (Diaz-Pulido and Garzón-Ferreira, 2002). Our findings are supported by the study of Meija et al. (2012), suggesting that macroalgae in reef environments typically show strong responses to fluctuations in water temperature and nutrient availability. Furthermore, Diaz-Pulido and Garzón-Ferreira (2002) detected similar variations in algal assemblages at the sheltered site of Chengue Bay and stated that this was mainly correlated to upwelling-related changes in water temperature although effects of nutrient availability could not be excluded. Since inorganic nutrient concentrations during upwelling events were generally above the reported thresholds, which sustained macroalgal blooms (dissolved inorganic nitrogen ∼1 µmol.L −1 ; phosphates ∼0.1 µmol.L −1 ; Lapointe, 1997), we suggest that nutrient availability is one of the main factors causing pronounced shifts in local algal assemblages.
Our results further showed that local bay-specific factors such as surface area or anthropogenic impact ( Table 1) could explain similar fractions of the variance in benthic groups and coral assemblages (splits 4, 6, 7, 11, and 14; Figure 3B and splits 2 and 3; Figure 4B, respectively). Surprisingly, ISA revealed that the separation of benthic groups was mainly based on variation in general coral abundance. These findings therefore suggest that local bayspecific factors influence corals more than the other investigated groups. A reason for the large influence of bay-specific factors on coral assemblages could be the relative longevity of corals compared to other monitored benthic groups such as algae, as this may lead to a cumulative effect of local disturbance events (Jackson, 1991). Accordingly, terrestrial run-off and associated sedimentation (Antonius, 1972), landslides (personal observations), dynamite fishing (Garzón-Ferreira and Cano, 1991) and other local disturbances in the past could have caused bay-specific differences in general coral abundance as well as coral species composition. In addition, the supply-side ecology of hermatypic corals may also explain the bay-specific differences in coral assemblages, as the variation of larval input, in general, determines the structure of local adult populations (Hughes et al., 2000). The relatively low dispersal ability of corals may lead to self recruitment and local differences of coral communities (Harrison and Wallace, 1990). For the study area, the mostly local effect of supply-side ecology on structuring coral assemblages is further enhanced by the fact that most corals in TNNP spawn when water surface currents are generally weak (i.e., August-October; Van Veghel, 1993;Acosta and Zea, 1997;de Graaf et al., 1999;Bayraktarov et al., 2014a).

LIMITATIONS OF THIS STUDY AND OUTLOOK
Our analyses indicated that wave exposure is the major factor structuring benthic communities in TNNP. However, given the calculated wave exposure for the investigated sites, hydrodynamic regimes differ from bay to bay as suggested by Garzón-Ferreira and Cano (1991) and exhibit site-specific temporal variations as shown by Bayraktarov et al. (2014b). Additional community assessments along the hydrodynamic gradient within the bays of TNNP could elucidate the effects of hydrodynamics on local benthic communities as proposed for coral assemblages by Werding and Sánchez (1989). One interesting aspect of the current study is that, although the overall coral abundances did not contribute to the differences in benthic group community structures between exposed and sheltered sites as depicted by the ISA analysis, species-specific differences in coral assemblages were present between sheltered and exposed sites. This raises the question about the power of using broad benthic groups (i.e., higher taxa) for revealing patterns in community structures. In fact, Diaz-Pulido and Garzón-Ferreira (2002) also found species-specific differences in algal assemblages at two study sites in the sheltered part of one TNNP bay, thus further emphasizing the need for detailed studies on lower taxonomic levels.

CONCLUSIONS
The current findings largely support our working hypothesis that wave exposure as regionally-controlled local parameter is the main driver of benthic community structure in the bays of TNNP. On the one hand, it causes specific community structures within exposed and sheltered sites across bays and, on the other hand, accounts for significant differences between exposed and sheltered sites within bays. Interestingly, upwelling events as important regional processes seem to play a less important role than wave exposure, particularly for coral assemblage structures. One of the reasons could be that seasonal upwelling events do not cause major spatial differences in water parameters in TNNP and mainly act at a temporal scale.
In summary, benthic community composition in TNNP is primarily driven by wave exposure (regionally-driven local factor), followed by upwelling (non-local factor) for major benthic groups, and bay-specific differences (local factors) for both benthic groups and coral assemblages.
The study not only highlights the existence of complex hierarchical levels of environmental factors acting on coral dominated communities, it also stresses the role of both local and non-local factors for understanding extant patterns in these highly complex ecosystems. We suggest that factors operating on different spatial scales should be integrated in future studies assessing community composition and resilience capacities of coral reefs.