Ecological patterns of macrofauna in sandy beaches of Costa Rica: A Pacific-Caribbean comparison

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

Contents lists available at ScienceDirect

Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss

Ecological patterns of macrofauna in sandy beaches of Costa Rica: A PacificCaribbean comparison

T

Jeffrey A. Sibaja-Corderoa,b,∗, Yolanda E. Camacho-Garcíaa,b, Juan Carlos Azofeifa-Solanoa, Bárbara Alvado-Arranzc a

Centro de Investigación en Ciencias del Mar y Limnología (CIMAR), Universidad de Costa Rica, Costa Rica Escuela de Biología, Universidad de Costa Rica, Costa Rica c Ciencias del Mar, Universidad Católica de Valencia - San Vicente Mártir, Spain b

A R T I C LE I N FO

A B S T R A C T

Keywords: Community ecology Coastal zone Environmental factors Intertidal invertebrates Tropical beaches

The present study compares the ecological patterns of macrofauna in sandy beaches between the Caribbean and Pacific coasts of Costa Rica. At each beach, the intertidal zone was divided into five strata from low to high tide level. In each stratum, sediment samples (corer area: 20.2 cm2, diameter: 5.07 cm) were collected to analyze the macrofauna. Additionally, grain size composition, total organic matter and carbonates of the sediments were determined. The macrofauna diversity was higher in the Pacific coast (15 taxa, per beach) than in the Caribbean (4 taxa, per beach). The most diverse taxon group was the annelids, followed by arthropods and mollusks. Within the mollusks, the bivalves and the gastropod family Olivellidae were only present in the Pacific coast. Both Caribbean and Pacific coasts had a vertical zonation of taxa distribution from high to low tide level. Typically, the low tide level was populated by several polychaetes, mid littoral by crustacean and mollusks, and the high tidal level was populated by isopods (Cirolanidae). Faunal differences in richness, abundance and composition of assemblages between both coasts can be explained by environmental differences. The tidal range is narrower on the Caribbean (0.5–1.5 m) than the Pacific coast (2–3 m). The slope of the beach is steep in several beaches of the Caribbean, whereas the Pacific has several dissipative beaches with gentle (or flat) slopes. The organic matter was lower (1.00%) in the Caribbean than in the Pacific (2.25%). In the Caribbean coast, the sediment is ∼90% fine sand (500–63 μm), while in the Pacific this fraction was ∼65%, resulting in more heterogeneous sediments in the Pacific. Finally, evidence of change in the abundance and species composition with the latitude in each coast was detected; indicating that these communities are highly variable within a local scale, which could be due to the diversity of sandy beaches and environmental variation that Costa Rica has in both coasts. Our results indicate that beaches with higher number of species or abundance were not necessarily designated as Marine Protected Areas. This information should be considered in the future for the establishment of new Marine Protected Areas.

1. Introduction Sandy beaches represent a multivariate space that is exposed to a dynamic physicochemical environment (Defeo et al., 2009). The sediment composition can vary within and between sites, resulting in dissimilar suitability to contain biological populations. Sandy beaches have a distinct fauna, well adapted to the harsh conditions of these habitats, such as desiccation and high temperatures during low tides, and the sediment instability and abrasion produced by waves during high tide and storms (Alongi, 1990). The biological assemblages at sandy beaches can change in different spatial and temporal scales (Defeo et al., 2009). The identity of species ∗

differs between and within broad geographical regions (Defeo et al., 2017). In a small scale, the highest variation occurs within the intertidal zone, both vertically and horizontally, resulting in a patchy distribution of biota. The species composition is influenced by the geographical and ecological gradients that determine the kind of animals that inhabit the sediments, and the biological interactions that regulates the population dynamics (Checon et al., 2018; Sibaja-Cordero, 2018). Barboza and Defeo (2015) and Defeo et al. (2017) presented the variation patterns of macrofauna richness for several sandy beaches around the world. However, most of the knowledge regarding the ecological patterns in sandy beaches refers to temperate and subtropical regions, rather than lower latitudes (< 10° N or S) (Defeo and

Corresponding author. Escuela de Biología, Universidad de Costa Rica, Costa Rica. E-mail address: jeff[email protected] (J.A. Sibaja-Cordero).

https://doi.org/10.1016/j.ecss.2019.04.032 Received 1 October 2018; Received in revised form 12 April 2019; Accepted 17 April 2019 Available online 25 April 2019 0272-7714/ © 2019 Elsevier Ltd. All rights reserved.

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

influence can expand the width of these beach zones in both coasts. The samples were collected in 12 and 16 sandy beaches of the Caribbean and Pacific coasts, respectively. The latitude and longitude of each beach were recorded. At each beach, five strata were established from low to high tide levels (I: Infralittoral, II: Low littoral, III: Mid littoral, IV: High Littoral and V: Sublittoral) as in Dexter (1974). At each stratum, five corers (20.2 cm2 area or 5.07 cm diameter, and 15 cm deep into the sediment) were introduced in the sediment to collect benthic fauna, distanced every 2 m. Each sample was preserved in 95% ethanol, and stained with Rose Bengal. In the laboratory, the samples were washed through a mesh sieve to retain the macrofauna (> 500 μm). The organisms were sorted from the sediment using a stereoscope (30x) and placed in vials with 95% alcohol. The fauna was classified to the most specific taxonomic category possible using available keys and guides for mollusks (Keen, 1971; Coan and Valentich-Scott, 2012); polychaetes (de León-González et al., 2009; Fauchald, 1977); decapod crustaceans (Fischer et al., 1995); and isopods (Brusca and Iverson, 1985). Voucher specimens were deposited at the Museo de Zoología, Universidad de Costa Rica. An additional sediment sample was obtained at each of the five strata during the sampling in all the beaches to measure the sediment variables. These sediment samples were dried at 60 °C using an oven. A subsample of 100 g was used for a granulometric analysis to determine the contribution of the different grain sizes (> 4000, 2000, 1000, 800, 710, 500, 300, 250, 125, 63 and < 63 μm). The total organic matter (TOM) in the sediments was determined using the loss on ignition method. Three 5 g subsamples were weighed in a ceramic crucible and burned in a stove during 4 h at 450 °C and weighed again for determination of TOM (Holme and McIntyre, 2005). The mean of the three subsamples correspond to the TOM percentage. The percentage of carbonates in the sediment was measured using one 5 g sediment subsample, which was dried at 90 °C in a previously weighed ceramic crucible for 1 h. After cooling in a desiccant jar, the sample was weighed. The sample and crucible were then placed in a muffle (550 °C for 1 h). After cooling, the sample was again weighed. Finally, the sample was returned to the muffle (1000 °C for 1 h). The weight loss between the 550 °C and 1000 °C represents the amount of CO2 lost from carbonate minerals, and the percentage of CaCO3 initially in the sediments was finally estimated (Dean, 1974). Other variables measured at each beach were the intertidal width (m), the mean slope (using a clinometer every 3 m), water salinity (PSU), and temperature (°C) of sand during the morning low tides. The beaches were classified regarding their conservation status: Marine Protected Area (MPA) or not; according to three categories of urbanization: urbanized (in the core of the town: > 50% of urbanized area), semi-urbanized (in residential areas: < 50% and > 30% of urbanized area) or natural (outside the core of the town: < 30% of urbanized area) (Villacampa et al., 2017).

McLachlan, 2013; Barboza and Defeo, 2015). The Central America region is a gap on the current knowledge about ecology of sandy beaches. Due to its isthmus condition, Central America represents an ideal region to analyze ecological data from sandy beaches at two different water masses (the Pacific Ocean and the Caribbean Sea). The results from this region can be used to understand the ecology of tropical sandy beaches, and to develop management measures in tropical beaches. An adequate understanding of the biological communities composition and distribution at tropical beaches will help to test the “habitat tropicalization” hypothesis, driven by climate change, using a worldwide distributed environment, such as the sandy beach (Barboza and Defeo, 2015). In Central America, Dexter (1972, 1974) studied several beaches on both Pacific and Caribbean coasts of Costa Rica and Panama, presenting data on the vertical zonation of invertebrates in the intertidal zone. The ecological studies of Dexter (1972, 1974) focused on richness, density and autoecology of the most common species. Other studies of beaches in Costa Rica have focused on biodiversity lists: Dittmann and Vargas (2001) presented a diagram of the faunal zonation in a beach within the Gulf of Nicoya; Fowler (1979) reported on insects and crabs in turtle nests at Tortugero beach; Sibaja-Cordero et al. (2014) and CorralesUgalde and Sibaja-Cordero (2015) showed the number of species in sandy beaches at the North and South Pacific coast, respectively; and Sibaja-Cordero (2018) presented data of spatial distribution of four species in the low littoral of a sandy beach in the North Caribbean of Costa Rica. In the studies mentioned above, environmental data is only presented as descriptive or qualitative information of some of these beaches. Quantitative ecology analyzing geographical distribution of biodiversity and macrofaunal assemblages and the environmental influence is necessary for tropical beaches. Beyond the studies of Dexter (1972, 1974), multivariate analyzes can help to understand the environmental factors driven the faunal diversity and composition in this habitat. This kind of scientific progresses is necessary for a better management of tropical sandy beaches. The present study compares the ecological patterns in sandy beaches between the Caribbean and Pacific coasts of Costa Rica. This study aims to (1) determine the differences of sediment and beach characteristics between both coasts and tidal levels; (2) analyze the macrofauna species richness and abundance considering both coasts, tidal levels, and environmental covariates; (3) determine the influence of environmental (sediment and beach characteristics, conservation and urbanization status) and spatial factors (coasts, beaches, and tidal levels) on explaining the similarity of macrofaunal assemblages. 2. Materials and methods 2.1. Methodology This study was carried out on sandy beaches of Costa Rica, Central America (Fig. 1). The area is subject to tropical climate, with ocean water temperature ranging usually between 22 and 30 °C. The Costa Rican Pacific coast (1254 km length) presents a dry season (Dec–Apr) and rainy season (May–Nov), however in the North Pacific area occurs a coastal upwelling during the dry season, driven by the NE winds, when water temperature drops below 20 °C (Fischer, 1981; Cortés, 2014). In the Pacific, the sandy beaches are interrupted by long sections of cliffs or rocky platforms. Several of the beaches are located within gulfs, bays, and coves, while others are exposed directly to heavy wave influence. The Caribbean coast has two rainy seasons through the year (Nov–Mar and Jun–Aug) (Cortés et al., 2010). The Costa Rican Caribbean coast is a straight line of 212 km, mostly composed of sandy beaches, interrupted by small recent reefs or emerged reefs by tectonics, mainly in the south Caribbean (Cortés et al., 2010). The tidal range is about 2.5–3 m in the Pacific, and at maximum 0.70 m (with a mean of 0.21 m) in the Caribbean coast (Lizano, 2006). Moreover, the wave

2.2. Environmental data To determinate the differences of environmental characteristics between the coasts and tidal levels, sediment data were analyzed (to obtain the asymmetry, selection, kurtosis, mean grain size, and the percentages of gravel, sand and fine sand) using the package ‘rysgran’ in R (Gilbert et al., 2012), following the method of Folk and Ward (1957). A Principal Component Analysis (PCA) biplot (Axis 1 and 2) was applied to the sediment data matrix, using the package ‘vegan’ in R. The other significant axes (eigenvalue > 1) are shown in supplementary material (Fig. A1). The goal of this method was to identify the highest correlated sediment fractions to reduce the data matrix that will be employed in the posterior analysis. The angle of separation between vectors of the variables in the biplot is proportional to their correlation (Legendre and Gallagher, 2001); this information and the loadings of the PCA axes (Table A1) were used to make groups of sediment fractions. Additionally, a ternary plot (supplementary material, Fig. A2) 95

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Fig. 1. Study site: Geographic position of A) Central America isthmus, and B) Costa Rica. C) Sandy Beaches in the Pacific Coast from North to South (blue circles): Ju: Junquillal, B: Playa Blanca, LC: Los Cocos, Eg: Estero Grande, Nn: Naranjo, N, Ns: Naranjo, S, Pe: Penca, Po: Potrero, Con: Conchal, PM: Punta Morales, Co: Cocorocas, P: Puntarenas, Cal: Caldera, J: Jacó, G: Golfito, Z: Zancudo, and in the Caribbean coast from North to South (green circles): W: Westfalia, V: Vizcaya, Bn: Bananito, N, Bs: Bananito, S, Ca: Cahuita, PV: Puerto Viejo, Pi: Pirripli, Co: Cocles, Ch: Chiquita PU: Punta Uva, M: Manzanillo, G: Gandoca. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Based in the principle of parsimony, the model was simplified using the step command in R (Crawley, 2007). The mean values ± 95% confidence of taxa and abundance by coast is presented.

was performed with the groups of sediment fractions that better discriminated the samples between beaches (results of the PCA), using the package ‘ggtern’ in R, to obtain groups of beaches with similar sediment composition. An Analysis of Variance by linear mixed-effects models, using the package ‘lmer’ in R (Crawley, 2007), was applied to each of the sediment fractions groups and sediment characteristics to evaluate differences of the mean value of each variable between coasts (Caribbean or Pacific), tidal level, and their interaction. The beaches were introduced in the model as a random effect. The same test was carried out with the TOM and carbonates data. For beach slope, salinity, and width of the intertidal; a two sample t-test was used to compare the mean value between coasts. All assumptions of the analysis were checked using the Shapiro-Wilk test for normality and the Bartlett test for homoscedasticity. The generalized variance inflation factor GVIF was calculated to check colinearity; a value of GVIF^(1/(2*Degres of fredom)))^2 less than 4 indicating no colinearity (Fox and Monette, 1992; Fox, 2002; O’Brien, 2007). The data was transformed by Box-Cox if necessary, using the package ‘MASS’ in R (Crawley, 2007). The percentage of the sediment fractions was expressed in radians for the analysis (Krebs, 1999).

2.4. Assemblages and environment A Hellinger transformation was applied to the abundance of each species (Legendre and Gallagher, 2001) because the data matrix presented many zeros, which represent a bias when the main goal is to study long spatial gradients (Legendre and Gallagher, 2001; Pec et al., 2017; Powell, 2018). Then, multivariate analysis (Euclidean-based methods) were applied to analyze the biological data matrix. All analyzes for assemblages were carried out in package ‘vegan’ in R. The transformed biological matrix was analyzed with a Multivariate Analysis of Variance Using Distance Matrices (adonis) to test if the assemblage composition changes between coasts and tidal levels, and their interaction (Powell, 2018). The beach identity was introduced in the analysis as a blocking factor to control the permutations of samples within each site (hierarchical spatial structure). The composition of assemblages was plotted using a PCA based on the transformed abundance data (Supplementary material, Fig. A3). Additionally, a SIMPER test based on euclidean distances of the transformed data was used to find the species that contributed the most to the change between coasts and between tidal levels within coasts. This analysis was done using the PAST software (Hammer et al., 2001). Variation Partitioning based on redundancy analysis (RDA) was used to measure the amount of variance of the transformed biological data that is explained by intertidal characteristics, sediment data, conservation and urbanization status, and spatial structure obtained by Principal Coordinates of Neighborhood Matrix (PCNM). The PCNM procedure was carried out to evaluate the difference in biological assemblages by different spatial scales independently of the environmental factors (e.g. Borcard et al., 2004). In some case this spatial patterns are result of the biology of the species as dispersal or competition process, and clumped or random distribution patterns (Borcard et al., 2004; Dray et al., 2006; Rodil et al., 2018). In this study, the PCNM resulted in 5 linear axes that can be used as spatial variables to correlate with biological matrix. First axes tend to represent large-scale patterns, while the later axes tend to represent smaller scale patterns

2.3. Abundance, richness, and environment To determinate the effect of the coast (Caribbean or Pacific), MPA status and the covariates latitude, urbanization degree, TOM, carbonates, sediment composition (mean values per beach), and their interaction with the coast factor over the species richness and total abundance, a generalized linear model was carried out (using the package ‘glm’ in R). The model for species richness was fitted using the Poisson family distribution, and for the abundance with Gamma family distribution, both with the log-link function, after an inspection of the fitted values, residuals and residual deviance (Crawley, 2007). Previously, to avoid singularity of covariates in the model, a Pearson index was used to identify high correlated variables (> 0.65). The most correlated covariates were: the sorting with TOM; the mean grain size with the gravel and fine sand; the coarse sand with the fine sand; and the MPA status with the urbanization degree (Table A2). From the previous list, only sorting, mean grain size, coarse sand, and MPA status were introduced in the model with the other uncorrelated variables. 96

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

(23–36%) in the Pacific, and an opposite decreasing pattern (39–18%) in the Caribbean (Tables 1 and 2). The fine sand percent increase from low to high tide level in the Caribbean, whereas in the Pacific the lowest values were at the mid tide level (Tables 1 and 2). Additionally, the contribution of coarse/sand is a good discriminator between beaches with low gravel contribution (< 10%) (Fig. A2). The sediment samples in the Pacific showed a trend towards poorly sorting (mean = 2.4, Fig. 2A), whereas the Caribbean had well sorted sediments (mean = 1.4, Tables 1 and 2). The fine sand samples result in more compact sediments (Table 1). The sorting is better in the high tide than in the low tide level in both coasts (Tables 1 and 2). Moreover, a difference between coasts was not detected with the mean values of kurtosis (Table 2). The mean value of TOM was higher (2.25%) in the Pacific than in the Caribbean (1.00%), and in the two coasts the value increased (∼0.65%) from high to low tide level (Tables 1 and 2). Both coasts presented a similar carbonates percentage mean (Tables 1 and 2). The higher values were found in the mid-levels of the intertidal (Tables 1 and 2). The temperature during low tide to high tide was of 31 °C and 36 °C, respectively, in the Pacific, and 28 °C and 33 °C in the Caribbean. The Caribbean coast had lowest water salinity (28.17 ± 1.27 PSU) than the Pacific (31.61 ± 1.63) (t = −3.26, df = 25.79, p = 0.003). The intertidal width was shorter in the Caribbean (18 ± 8 m, ± 95% CI) than in the Pacific (99 ± 60 m)(log(x), t = −5.61, df = 20.7, p < 0.001). This difference resulted in a highest mean slope in the Caribbean (5.88 ± 1.49%) than in the Pacific (3.62 ± 1.15%) (t = 2.35, df = 22.15, p = 0.028).

Fig. 2. PCA (scaling type 2) of correlation between fractions and statistics of the sediment samples (eigenvalues: PC-1 = 5.8660 and PC-2 = 3.1225), from the Pacific and Caribbean coast of Costa Rica. PC-3 = 2.0130 and PC-4 = 1.5707 are shows in Fig. A1.

(Pec et al., 2017; Powell, 2018). With the present sampling design, the PNCM objects represent distances between and within coasts. PNCM-1 axis indicated a separation between the Caribbean and Pacific beaches, and between the North and South Pacific (Fig. A4); PNCM-2 indicated the separation of Jacó Beach; PNCM-3 indicates that the beaches in Gulf of Nicoya and Golfo Dulce (also with PNCM-4) are different spatial regions than the other beaches. In PNCM-5 the presence of gradients from northern to southern beaches in North Pacific and Caribbean coast are suggested (Fig. A4). Each partition obtained in Variation Partitioning was tested with partial RDA using an ANOVA test, and plot in a Venn diagram. The axes and general model of RDA were tested with an ANOVA, and variables in the model were checked for collinearity with a VIF (value < 10 is the threshold in Table A4) (Pec et al., 2017; Powell, 2018).

3.2. Abundance, richness, and environment A total of 1266 individuals of 191 taxa distributed in 12 phylla were found in the 28 sandy beaches sampled in the present study (Table A3). The arthropods (mainly crustaceans) accounted for 42.5% of the total abundance and were represented by 47 taxa, followed by polychaetes with 41.3% of abundance in 89 taxa, and mollusks with 5% of abundance with 26 taxa (Table A3). The mean number of taxa per beach was higher in the Pacific coast (15 ± 9 taxa per beach with a maximum of 68) than in the Caribbean (4 ± 1 taxa per beach with a maximum of 6). This difference can be explained considering the interaction of the coast with several covariates, such as the negative relationship of taxa richness with the beach slope, positive with the intertidal width (Fig. 3), and slightly negative with the carbonates content (Table 3, Fig. 3). Otherwise, positive relationships were found between taxa richness and the mean grain size, and the kurtosis of the sediments (and indirectly with the sorting, TOM), that changed the value of their pendants between coasts (Table 3, Fig. 3). The taxa richness decreased with the positive skewness (to the fine sediments). The MPA status shows higher taxa values in the Pacific coast than in the Caribbean, but the highest values were found in non MPA sites (Table 3, Fig. 3). Variables that had no effect on the species richness were removed from the model (Table 3). The total abundance was highly variable between coasts, with a density of 635 ± 472 ind. m−2 in the Pacific and 1090 ± 551 ind. m−2 in the Caribbean (Table 4). In the Caribbean the slope was above 6°, with a narrow width in the intertidal of several beaches and the abundance reached lower values at these sites (Table 4, Fig. 4). The abundance was highly variable in the Pacific, but presented high abundance in beaches with the widest intertidal available. Additionally, the abundance decreased in the Caribbean when the sediments had higher mean grain size and coarse sand contribution, but increased with the grain size in the Pacific (Table 4, Fig. 4), which resulted in a decrease of abundance with positive (fine) skewness in the Pacific coast, but not in the Caribbean (Table 4, Fig. 4). The abundance increased with the sorting, kurtosis and carbonates in the Pacific and decreased in the Caribbean (Table 4, Fig. 4). Finally, the abundance showed a decrease with the increase of salinity (Table 4, Fig. 4); and a change of the

3. Results 3.1. Environmental data The sediment data resulted in three groups of high correlated fractions in the PCA (Fig. 2, Fig. A1, Table A1). Gravel sand comprised all the fractions > 1000 μm, the coarse sand with fractions < 1000 μm and > 250 μm, and fine sand (< 250 μm). Additionally, the median and mean grain sizes were highly correlated with the sorting, and with the gravel sand (Fig. 2). The skewness and kurtosis of sediments were positively correlated between each other and with the 250 μm fraction. These variables were negatively correlated with the mean grain size and gravel sand (Fig. 2). The gravel, the coarse, and fine sand produced the best separation of the samples from the beaches of Costa Rica. In the ternary diagram of supplementary material, the sediments of the Pacific coast of Costa Rica presented high contribution of grain size > 500 μm than in the Caribbean (Fig. A2), but the beaches within each coast presented high variation in their sediments (Fig. A2). In general, Pacific beaches had higher grain size by higher mean gravel content compared to the Caribbean (21% v 1% in mean, respectively). This pattern occurred in all the tidal levels (Tables 1 and 2), producing a coarse skewness (negative values). The coarse sand content had similar general means between coasts (∼30%). However, this fraction showed an increasing value from low to high tide level 97

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Table 1 Mean and Confidence Interval (95%) of statistics and main fractions of the sediment samples by intertidal levels from the sandy beaches of the coasts of Costa Rica. Tidal Level Variable

Coast

General Mean

I: Low tide

II

III

IV

V: High tide

Mean grain size (μm)

Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific Caribbean Pacific

258 ± 27 543 ± 122 258 ± 27 641 ± 174 1.454 ± 0.051 2.43 ± 0.18 1.10 ± 14.00 21.60 ± 117.00 27.73 ± 51.93 30.08 ± 74.07 71.12 ± 85.93 48.32 ± 152.01 1.00 ± 0.01 1.07 ± 0.03 −0.014 ± 0.003 −0.018 ± 0.039 1.04 ± 0.08 2.25 ± 0.55 24.21 ± 5.70 18.13 ± 4.44

316 ± 76 538 ± 345 314 ± 74 691 ± 544 1.61 ± 0.20 2.54 ± 0.59 3.82 ± 3.60 20.85 ± 14.1 38.69 ± 17.75 22.71 ± 9.97 57.49 ± 19.40 56.44 ± 15.77 1.01 ± 0.07 1.09 ± 0.15 0.019 ± 0.064 −0.003 ± 0.190 1.44 ± 0.24 2.61 ± 0.71 31.59 ± 13.82 25.03 ± 13.09

264 ± 63 537 ± 301 262 ± 61 636 ± 487 1.45 ± 0.08 2.38 ± 0.44 0.98 ± 1.10 20.46 ± 13.80 30.33 ± 19.77 27.26 ± 10.27 68.72 ± 20.45 52.28 ± 15.12 0.99 ± 0.09 1.16 ± 0.13 0.004 ± 0.056 0.074 ± 0.188 1.02 ± 0.28 2.28 ± 0.72 26.02 ± 15.04 20.31 ± 11.29

238 ± 48 618 ± 291 238 ± 47 677 ± 361 1.42 ± 0.05 2.42 ± 0.48 0.33 ± 0.26 25.44 ± 14.27 25.57 ± 17.88 33.83 ± 10.43 74.1 ± 18.08 40.73 ± 13.62 1.05 ± 0.10 1.06 ± 0.18 0.024 ± 0.050 −0.006 ± 0.19 1.04 ± 0.22 2.13 ± 0.52 12.50 ± 12.55 19.49 ± 10.00

251 ± 62 577 ± 216 251 ± 61 706 ± 299 1.39 ± 0.04 2.61 ± 0.55 0.50 ± 0.53 25.66 ± 12.34 25.67 ± 19.62 36.21 ± 8.36 73.83 ± 20.12 38.13 ± 11.91 0.95 ± 0.08 1.02 ± 0.16 −0.030 ± 0.049 −0.100 ± 0.168 0.93 ± 0.17 2.3 ± 0.73 18.59 ± 11.16 13.42 ± 7.68

223 ± 47 445 ± 181 225 ± 47 496 ± 235 1.41 ± 0.05 2.20 ± 0.42 0.13 ± 0.16 15.60 ± 11.35 18.41 ± 15.80 30.38 ± 8.38 81.46 ± 15.92 54.02 ± 13.79 1.02 ± 0.07 1.02 ± 0.12 −0.09 ± 0.048 −0.05 ± 0.143 0.79 ± 0.23 1.93 ± 0.57 22.20 ± 11.82 12.39 ± 5.86

Median grain size (μm) Sorting (Geometric scale) Gravel (> 1000 μm) (%) Coarse sand (< 1000 - > 300 μm) (%) Fine sand (< 300 μm) (%) Kurtosis (Geometric scale) Skewness (Geometric scale) TOM (%) CaCO3 (%)

related to the higher grain size and sorting. Quadrant A and D in RDA contains most of the northern beaches in both coasts, and in quadrant B and C are mainly beaches in southern direction of both coasts (Fig. 5). The macrofauna assemblages were associated with an increase in carbonates in South Caribbean beaches (quadrant B and C in RDA), and with the increase in gravel, TOM and skewness in the Pacific beaches (quadrant B in RDA). The beaches within the estuarine zones of Gulf of Nicoya and Golfo Dulce with wide intertidal and higher TOM values are clustered in quadrant B of the RDA (Fig. 5). A trend for a vertical zonation within the intertidal was also detected using the RDA. In quadrant B, 29% of the samples are from level I, 25% from level II (Fig. 5). Likewise, in quadrant C, 25% of the samples are from level I, and 33% from level II. In quadrant D of RDA a 47% of the 17 samples comes from level II and III (mid littoral). Finally, quadrant A contains several samples from the higher intertidal zone, specifically from level III (26%), IV (39%), and V (30%). Moreover, evidence of vertical zonation was not found for seven of the 26 beaches by the higher similarity between samples. Six of these beaches were from the Pacific coast and are placed in the quadrant B of the RDA ordination (Fig. 5). The Variation Partitioning in the Venn diagram (Fig. 6) indicates that only 19% of the variation of macrofaunal assemblages can be explained by the factors taken into account in the present study. The geographical position of the beach and the sediment characteristics were the main drivers influencing the composition of animals in the sandy beaches, this due to the fact that each matrix explains > 0.04 of the total variation with an r2 = 0.10 (Fig. 6). The beach features, such as intertidal width, beach slope, and tidal level, occupied the third place explaining the biological assemblages. The MPA and urbanization status per se did not explain the macrofaunal assemblages (partial RDA, variance = 0.01, residual variance = 0.48, F = 1.24, p = 0.177), and only increased their influence in relation with geographic location, sediment type, and beach profile (Fig. 6). The partial RDA detected an effect of the partitions showed in the Venn diagram over the faunal composition (Table A5). These matrices acted in synergy to explain the synecology of the macrofauna (r2 = 0.19).

abundance occurs with the latitude (Table 4, Fig. 4). Variables that had no effect on the abundance were removed from the model (Table 4). 3.3. Assemblages and environment The macrofaunal assemblages varied in their composition between Caribbean and Pacific coasts (Table 5, Fig. A3). The main difference between coasts was accounted by the high abundances of the isopod of the species complex Excirolana braziliensis, the polychaetes Scolelepis squamata, and Pisionidens indica in the Caribbean in comparison to the Pacific coast (Table 5). A high number of rare species were observed in the Pacific (Table A3). Mollusks (gastropods and bivalves) were only found in the Pacific coast and contributed with 5% of dissimilarity between coasts (Table 6). The observed macrofauna follows a vertical zonation pattern within the intertidal zone at each coast (Table 5, Fig. A3). This pattern, however, was not similar between coasts (Table 5). In the Caribbean, the infralittoral was dominated by the worm P. indica, followed by the spionid S. squamata in the lower littoral; and small nemerteans in the mid littoral (Table 6). The isopod E. braziliensis, the amphipods Talitridae and adult insects (Psocoptera) had more individuals in the high littoral and supralittoral levels (Table 6). These groups accounted for 85% of the differences between tidal levels. In the Pacific coast, the infralittoral had more abundance of the worm P. indica, the mole crab Emerita rathbunae, and several nematods. The polychaete Hemipodia simplex was more abundant in the low littoral, while in the mid littoral the worm Pisione remota and the gastropod Olivella semiestriata were present (Table 6). Finally, the isopod E. braziliensis, and adult insects of the group Staphylinidae had high abundances in the high littoral and supralittoral levels (Table 6). These species accounted for only 52% of the change between tidal levels. The difference in species identity resulted in highest dissimilarity of low and mid littoral zones between the coasts (II and III levels in Fig. A3). However, these spatial-scale factors explained only 13% of the variance. In the RDA ordination (Fig. 5), controlling the biological matrix by the spatial and environmental data (Model variance = 0.22, residual variance = 0.47, F = 2.70, d.f. = 21/118, p = 0.001, and VIF < 10, Table A4), also showed the separation of the macrofaunal assemblages between Pacific (quadrant A and B) and Caribbean beaches (quadrant C and D). The faunal composition in Caribbean beaches was positively related to the positive kurtosis by the increase of the fine sand fraction (Fig. 5). In the Pacific coast the macrofauna assemblages were more

4. Discussion In the present study, the diversity of beach fauna was higher in the Pacific than in the Caribbean coast. Moreover, several rare taxa occurred in low abundances in several of the sites along the Pacific coast. 98

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Table 2 Analysis of Variance using linear mixed-effects models with Satterthwaite's method of selected variables (based on PCA) of sediments samples from the sandy beaches of the coasts of Costa Rica. Gravel %

Sum. Sq

df

F

p

GVIF^(1/(2*Df))

Model: log (fraction in radians + 1) ∼ Coast * Level + (1 | Beach) Coast 0.136 0.136 Tidal Level 0.116 0.029 Coast*Tidal Level 0.088 0.022 Random effects: Beach Variance: 0.026

1/26 4/104 4/104 Residuals:

17.14 3.63 2.77 0.008

< 0.001 0.008 0.031

1.111 1.528 1.546

Coarse sand %

Sum. Sq

Mean. Sq

df

F

p

GVIF^(1/(2*df))

Coast Tidal Level Coast*Tidal Level Random effects: Beach

0.007 0.054 0.350 Variance:

0.007 0.014 0.088 0.031

1/26 4/104 4/104 Residuals:

0.88 1.66 10.78 0.008

0.357 0.163 < 0.001

1.095 1.528 1.543

Fine sand %

Sum. Sq

Mean. Sq

df

F

p

GVIF^(1/(2*Df))

Coast Tidal Level Coast*Tidal Level Random effects: Beach

0.153 0.452 0.900 Variance:

0.153 0.113 0.225 0.114

1/26 4/104 4/104 Residuals:

5.31 3.92 7.81 0.029

0.029 0.005 0.000

1.092 1.528 1.543

Mean. Sq

df

F

p

GVIF^(1/(2*Df))

0.200 0.013 0.016 0.009

1/26 4/104 4/104 Residuals:

41.73 2.69 3.38 0.005

< 0.001 0.035 0.012

1.178 1.528 1.557

Model: (1/Y2) ∼ Coast * Level + (1 | Beach) Sorting (Geometric scale) Sum. Sq Coast Tidal Level Coast*Tidal Level Random effects: Beach Model: fourth root (Y) ∼ Coast * Level

0.200 0.052 0.065 Variance: + (1 | Beach)

Mean. Sq

Kurtosis (Geometric scale)

Sum. Sq

Mean. Sq

df

F

p

GVIF^(1/(2*Df))

Coast Tidal Level Coast*Tidal Level Random effects: Beach

0.002 0.009 0.008 Variance:

0.002 0.002 0.002 0.001

1/26 4/104 4/104 Residuals:

0.65 0.92 0.87 0.002

0.428 0.452 0.485

1.536 1.528 1.618

TOM %

Sum. Sq

Mean. Sq

df

F

p

GVIF^(1/(2*Df))

Coast Tidal Level Coast*Tidal Level Random effects: Beach

0.152 0.240 0.028 Variance:

0.153 0.060 0.007 0.014

1/26 4/104 4/104 Residuals:

18.28 7.18 0.84 0.008

< 0.001 < 0.001 0.502

1.196 1.528 1.560

CaCO3%

Sum. Sq

Mean. Sq

df

F

p

GVIF^(1/(2*Df))

Coast Tidal Level Coast*Tidal Level Random effects: Beach

0.069 0.093 0.104 Variance:

0.069 0.523 0.026 0.218

1/26 4/104 4/104 Residuals:

0.53 4.06 0.20 0.129

0.712 0.004 0.936

1.193 1.528 1.560

reflective beaches (Barboza and Defeo, 2015). The tidal range is a good predictor of species diversity of sandy beaches (Barboza and Defeo, 2015; Defeo et al., 2017). This hypothesis is also supported by our results, where the tidal range is 2.5 m–3 m in the Pacific and only a maximum range of 0.7 m in the Caribbean (Lizano, 2006). Dexter (1972, 1974, 1976, 1979) found the same differentiation patterns, with higher richness in the Pacific than in Caribbean, in Mexico, Colombia, Panama, and Costa Rica. However, Defeo et al. (2017) found a higher diversity in the Atlantic beaches than in their Pacific counterpart in South America. This difference is possibly associated with the higher primary productivity in Atlantic beaches (Defeo et al., 2017). Ortega-Cisneros et al. (2017) show that input of organic matter from estuarine systems can influence the fauna in sandy beaches of South Africa. In Central America, the Pacific coast has higher productivity than the Caribbean, due to the presence of seasonal tropical upwellings (Tehuantepec, Papagayo, and Panama) and several estuarine systems that contribute with nutrient inputs (Dittmann and Vargas, 2001; Cortés et al., 2010). The Caribbean beaches had a narrower intertidal width than in the

Our results are in concordance with Checon et al. (2018), who found exclusive species on almost every sampled beach in Brazil. Moreover, in both studies (Checon et al., 2018; and present study), typical negative relationship between taxa richness and the beach slope was also observed (Defeo and McLachlan, 2013). This trend was evident in the Pacific of Costa Rica, where a wide range of beach slopes can be found. The range of reflective to dissipative beaches allows a diverse beach biota. In dissipative beaches, more gentle conditions and microhabitats occur in contrast with reflective beaches (Barboza and Defeo, 2015), by the habitat availability influence (Defeo et al., 2017). Otherwise, the negative taxa-beach slope association was not clearly detected in the Caribbean of Costa Rica, where most beaches had steeper slopes (mean 5.88°). Considering data from both coasts, this general pattern was overall, in agreement to Defeo and McLachlan (2013). In the present results, the dissipative beaches in the Pacific are one extreme, while the Caribbean reflective beaches represent the other one within this gradient. The few number of species living in the sediments of Southern Caribbean beaches are a non-random selection of organisms, which can inhabit these harsh environments due to their particular adaptations to 99

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Fig. 3. Response of number of taxa of macrofauna with the beach and sediments characteristics, of each sandy beach from Caribbean and Pacific coasts of Costa Rica. Broken line is the linear regression of taxa with the independent variable, blue and green line for the Pacific and the Caribbean coast, respectively. Similar regression slope between coasts: black broken line. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

scenario to scenario within a latitudinal zone. In the present study, the abundance and species composition of the fauna correlated with the latitude is the result of local environmental conditions. In Costa Rica, the beaches in the north Pacific are exposed to a seasonal upwelling, and the beaches in the gulfs (Nicoya and Golfo Dulce) are sheltered, have a low slope, heterogeneous sands and high nutrients (Dittmann and Vargas, 2001; Cortés, 2014). Sites between these systems are intermediate or reflective beaches with gravel and coarse sands. Environmental variables that changed between these zones of the Pacific coast are the variations in TOM and sediment grain size. In the Caribbean, the beaches change from wide intertidal with small slope degree, fine terrigenous sediments in the north, to narrow intertidal, steeper slope, high carbonate and coarse grains in the south (SibajaCordero, 2018; Cortés et al., 2010). Alongi (1989) mentioned the influence on the composition of fauna, by the reduction in food availability (TOM) in high carbonate sands. In this way, the higher variation in the composition of species was caused by the differences of the beach environment depending of their geographic position at each coast. The difference between coast assemblages is also reflected in the vertical zonation of organisms, because sediment distributions changed

Pacific, due to a different tidal range. Dexter (1974) considered this to explain the remarkable difference in abundance between coasts in Costa Rica. This difference was not found as a general rule in the present study, because the density of individuals resulted in higher variance. The populations of several species contributed for the total abundance in the Pacific coast, while some beach systems in the northern Caribbean coast can compensate the low diversity with fewer species that monopolized the total abundance. The low abundance was related to the narrowed width of the intertidal and the fine skewness and compaction of the carbonate in the southern Caribbean. Additionally, in the Pacific the increase in grain size, heterogeneity of sediment and intertidal width was associated to an increase in the abundance of macrofauna. Heterogeneous sediments (higher sorting value) contain more microhabitats for the organisms and beaches, and then can harbour higher densities of fauna (McArthur et al., 2010). Defeo and McLachlan (2013); Barboza and Defeo (2015); Defeo et al. (2017) found that latitudinal changes can explain the richness at regional and worldwide scales. Alongi (1990) explains that within a latitudinal zone, variations are greater than the differences among latitudes. According to this, the assemblages and their species vary from

Table 3 Simplified Generalized linear models for taxa number per sandy beach in relation to the environmental data by coast (Pacific and Caribbean of Costa Rica). Simplified model: Taxa number ∼ Coast + Intertidal width + MPA + Slope + CaCO3 + Mean grain size + Skewness + Kurtosis + Coast*Slope + Coast*CaCO3 + Coast*Mean + Coast*Kurtosis, family = poisson (log = link) Null deviance: 358.505 on 27 degrees of freedom Residual deviance: 11.904 on 15 degrees of freedom AIC: 138.26 Coefficients:

Estimate

Std. Error

z

p

Intercept Coast (Pacific vs. Caribbean) Intertidal width MPA Slope CaCO3 Mean grain size Skewness Kurtosis Coast*Slope Coast*CaCO3 Coast*Mean grain size Coast*Kurtosis

0.52 −2.17 0.004 −0.62 0.03 −0.01 0.005 −2.20 −0.23 −0.26 −0.02 0.003 4.70

1.56 2.03 0.001 0.25 0.10 0.01 0.001 0.54 1.64 0.12 0.02 0.002 1.91

0.33 −1.07 −5.85 −2.50 0.26 −0.48 2.67 −4.07 −0.14 −2.14 −1.43 −1.65 2.47

0.740 0.290 < 0.001* 0.010* 0.800 0.630 0.010* < 0.001* 0.890 0.030* 0.150 0.100 0.010*

*p < 0.05. 100

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Table 4 Simplified Generalized linear models for the total abundance per sandy beach in relation to the environmental data by coast (Pacific and Caribbean of Costa Rica). Simplified model: Abundance+1 ∼ Coast + Intertidal width + PSU + CaCO3 + Coarse +Mean + Sorting + Skewness + Kurtosis + Latitude + Coast*Intertidal width + Coast*PSU + Coast*CaCO3 + Coast*Coarse + Coast*Mean + Coast*Sorting + Coast*Skewness + Coast*Kurtosis, family = Gamma (link = log) Dispersion parameter: 0.3716762 Null deviance: 37.3370 on 27 degrees of freedom Residual deviance: 5.1806 on 9 degrees of freedom AIC: 249.28 Coefficients:

Estimate

Std. Error

z

p

(Intercept) Coast (Pacific vs. Caribbean) Intertidal width PSU CaCO3 Coarse Mean Sorting Skewness Kurtosis Latitude Coast*Intertidal width Coast*PSU Coast*CaCO3 Coast*Coarse Coast*Mean Coast*Sorting Coast*Skewness Coast*Kurtosis

14.37 −28.56 0.04 −0.46 0.10 −0.35 0.11 −15.59 −20.01 −5.12 0.74 −0.05 0.55 −0.13 0.33 −0.11 16.43 18.05 11.42

7.04 8.00 0.02 0.17 0.03 0.12 0.04 3.93 8.19 2.76 0.30 0.02 0.20 0.04 0.12 0.04 3.95 8.25 3.28

2.04 −3.57 2.04 −2.74 2.96 −3.00 2.85 −3.97 −2.44 −1.86 2.48 −2.43 2.76 −3.28 2.79 −2.82 4.16 2.19 3.48

0.071 0.006* 0.071 0.023* 0.016* 0.015* 0.019* 0.003* 0.037* 0.096 0.035* 0.038* 0.022* 0.010* 0.021* 0.020* 0.002* 0.056 0.007*

*p < 0.05.

vertically in different ways between both coasts. Celentano et al. (2019) found an increase in diversity and density from the dune zone to the swash zone of a dissipative beach, which was also observed in the present study. As mentioned above, high sorting values indicate high heterogeneous sediments, and in the present study the taxa richness and density increased in these kinds of sediments. The groups found in low littoral levels of both coasts were mainly small motile carnivore worms such a pisionids, glycerids, nepthyids, and nemerteans that can use the interstitial space of these heterogeneous sands, or they can burrow in the less compacted terrigenous sands using the labile material or feeding on other macrofaunal and meiofaunal species (Jumars et al., 2015). Other small forms found in the lower tidal levels of isolated beaches of the Pacific coast were Dorvilleidae, Protodrilidae, Polygordidae, and Saccocirridae. These organisms are also dependent of

Table 5 Multivariate Analysis of Variance using Distance Matrices of macrofauna from the sandy beaches by tidal level and coasts of Costa Rica.

Coast Level Coast:Level Residuals Total

Sum. Sq

Mean. Sq

df

F

p

r2

4.67 4.72 2.82 83.70 95.91

4.67 1.18 0.71 0.64

1 4 4 130 139

7.25 1.83 1.10

0.001 0.001 0.048

0.05 0.05 0.03 0.87 1

Blocks: Beach.

Fig. 4. Response of abundance of macrofauna with the slope of the beach and sediments characteristics, of each sandy beach from Caribbean and Pacific coasts of Costa Rica. Broken line is the linear regression of taxa with the independent variable, blue and green line for the Pacific and the Caribbean coast, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 101

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Table 6 SIMPER analysis to show the main macrofaunal taxa in contribution (Contrib. %) to the dissimilarity (Av. Dissim.) between sandy beaches of each coast (Caribbean and Pacific), and between intertidal levels (I = low tide, V = high tide) within each coast. The contribution of the main phyla is also presented. The density of intabledividuals m−2 ± 95% confidence is presented by coast and tidal level. Taxon

Av. Dissim

Contrib. %

Cum. %

Caribbean

Pacific

Phyllum

Total Contrib. %

Excirolana braziliensis Scolelepis squamata Pisionidens indica Hemipodia simplex Nemertea indet. I Talitridae indet.

0.36 0.26 0.07 0.05 0.04 0.04

25.39 18.59 5.048 3.855 2.926 2.926

25.39 43.99 49.03 52.89 55.82 58.74

328 ± 175 257 ± 158 13 ± 10 0 7±8 5±6

301 ± 264 5±6 15 ± 13 50 ± 30 0 0

Annelida Arthropoda Mollusca Nemertea Platyhelminthes Others

43 40 5.8 5.2 1.7 4.3

Caribbean

Av. Dissim

Contrib. %

Cum. %

I

II

III

IV

V

Excirolana braziliensis Scolelepis squamata Pisionidens indica Nemertea indet. 1 Talitridae indet. Psocoptera indet. A

0.39 0.35 0.10 0.08 0.08 0.03

32.01 28.90 8.575 6.704 6.589 2.746

32.01 60.90 69.48 76.18 82.77 85.52

322 ± 347 148 ± 223 33 ± 36 0 0 0

635 ± 624 610 ± 525 25 ± 25 8 ± 16 8 ± 16 0

214 ± 169 437 ± 488 8 ± 16 25 ± 35 0 0

380 ± 475 82 ± 114 0 0 8 ± 16 0

90 ± 69 8 ± 16 0 0 8 ± 16 8 ± 16

Pacific

Av. Dissim

Contrib. %

Cum. %

I

II

III

IV

V

Excirolana braziliensis Hemipodia simplex Pisione remota Nematoda spp. Olivella semiestriata Pisionidens indica Emerita rathbunae Nemertea indet. A Polycladida indet. C Staphylinidae indet. A

0.31 0.10 0.07 0.07 0.05 0.04 0.04 0.03 0.03 0.03

22.02 6.808 4.910 4.882 3.623 2.919 2.606 2.069 1.763 1.763

22.02 28.82 33.73 38.62 42.24 45.16 47.76 49.83 51.60 53.36

12 ± 17 68 ± 84 25 ± 38 105 ± 143 37 ± 72 50 ± 53 6 ± 12 31 ± 34 0 0

43 ± 85 74 ± 84 43 ± 85 105 ± 171 18 ± 36 6 ± 12 18 ± 36 31 ± 34 6 ± 12 0

217 ± 171 43 ± 43 49 ± 97 31 ± 49 18 ± 36 18 ± 36 0 25 ± 38 0 0

507 ± 663 62 ± 87 6 ± 12 167 ± 327 6 ± 12 0 0 0 0 0

730 ± 1121 0 0 0 0 0 0 0 0 6 ± 12

interstitial spaces and they cannot be found in the finest beaches of Costa Rica. Other finding of Celentano et al. (2019) reports that most organisms occurs in the surface layer of the sediments, and are small-sized. In the present study, the majority of the organisms found in the lower littoral have small bodies, and inhabit the superficial layer of the sediments. In the present study, the spionid S. squamata occupies the compacted carbonate sands of the lower littoral at Caribbean beaches. This is concordant with Alongi (1989), who explained that infaunal organisms are generally small surface deposit and suspension feeders in the tropical benthos. Macrofauna with larger body size was only found at the beaches close to mangrove systems in the Pacific coast. In such localities, the trophic webs are maintained by the constant input of organic matter (Dittmann and Vargas, 2001). These localities have also the most gravelliest sediments, and therefore relatively large organisms, such as fiddler and pea crabs (Ocypodids and Pinnotherids, respectively) were found at these beaches, as well as bivalves of the families Tellinidae, Donacidae, Ungulinidae, and Veneridae, occupying the low and mid littoral zones. The mole crab E. rathbunae (low littoral) and the cowrie O. semiestriata (mid littoral) were also conspicuous organisms, which were found in some clean fine sandy beaches of the Pacific coast. Both of these species have adaptations to bury into fine sediments and their presence in these beaches was expected (Troost et al., 2012; CorralesUgalde and Sibaja-Cordero, 2015). Most of the organisms mentioned above are not easily collected using cores, and future sampling with quadrants will be needed to better understand their densities. Small adults of insects inhabited the high littoral levels of the beaches (Psocoptera in the Caribbean and Staphylinidae in the Pacific). The presence of these groups is not a surprise in this habitat. Frank and Ahn (2011) recorded several adults of Staphylinidae from the upper intertidal of several sandy beaches around the world, and Colombini et al. (2005), found Psocoptera using traps in dunes of sandy beaches in

Morocco. The group of insects associated to sandy beaches is a topic that needs further research. Mourglia et al. (2015) found a relationship of certain families of Coleoptera with the finest sands in a Uruguayan beach. Sediments in high littoral levels of Costa Rica were finest in the Caribbean than in the Pacific, but trying to explain the presence of insects groups with this variable will need to be tested in the future. Other habitants of high tidal levels were isopods (E. braziliensis) and amphipods (Talitridae). The species complex E. braziliensis lives in all tidal levels but adults are frequently found in the high tide level (Glynn et al., 1975). These organisms bury in the sand during the low tide, as a strategy to tolerate high temperatures and desiccation, while during high tide they come out of their burrows and feed on other invertebrates and carrion (Brusca and Iverson, 1985). Checon et al. (2018) found a variation in taxa composition within the same coast in Brazil, and Rodil et al. (2018) found that composition and abundance changed in the North coast of Spain by environmental and spatial scale factors. As in the present study, the composition of organisms was mainly related to spatial variables in first place, followed by sediment characteristics, and beach profile. Although the previous factors explained the composition of beach fauna, there was also a high percentage of variance that was not explained by the models. The beach fauna has different strategies for their spatial distribution at different scales (Sibaja-Cordero, 2018; Rodil et al., 2018). Differential recruitment, intraspecific and interspecific competition are some factors that can explain this variation (Borcard et al., 2004; Dray et al., 2006; Checon et al., 2018). The results are consistent with Defeo and McLachlan (2013) for a macroscale comparison between beaches. The authors point out that more sampling is needed in sandy beaches in order to understand the species-area relationship. The results of the present study indicate that the MPA and urbanization degree do not play a relevant role to increase the number of species. The MPAs in the Caribbean coast presented lower taxa number and densities than beaches without protection status, because the 102

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

dissipative beaches are in the non-protect areas along this coast. Moreover, in the Pacific coast the numbers of species did not differ between MPA status or urbanization degree. This is similar to previous assessments of sandy beaches in the Northern Pacific of Costa Rica, because beaches with higher number of species or abundance were not necessary designated as MPAs (Sibaja-Cordero et al., 2014). During the present study beaches were sampled once, however is necessary, to develop a monitoring plan for sandy beaches in both coasts to test the effect of MPAs on the management of sandy beaches biodiversity. Nonetheless, is clear that several MPAs have been established to protect others habitats, manage fisheries, or protect recreational areas. Considering the present analysis based on the composition of the fauna, the MPA status would play a relevant role when its design, establishment, and management are taken into account along with the geographical area, sediment type and beach profile. For example, in the case of Northern Pacific MPAs (most dissipative beaches), higher values of density and taxa were found in comparison with the Caribbean MPAs (most reflective beaches). Defeo et al. (2009) points out the critical need of taking into account a more comprehensive ecological understanding of sandy beaches when designating coastal MPAs.

5. Conclusions The number of species found in this study was lower in the Caribbean than in the Pacific coast, related to the narrower tidal range, higher sediment compaction and steeper beach slope in the Caribbean beaches than in the Pacific. Similar to previous studies in this habitat, the spatial factors and sediment characteristics were main drivers of the macrofauna abundance and composition. Additionally, the fauna in sandy beaches was correlated with latitude, because of the interaction of geographic and local environment influencing the sediments composition and beach profile. Further studies on the biological interactions are needed to understand the high variation not associated to the analyzed factors in this study. The establishment of MPAs including sandy beaches and shorelines needs to consider these ecological aspects when aiming to preserve the species diversity and the adequate conservation and management of these coastal habitats.

Fig. 5. Redundancy analysis (RDA) ordinations (scaling type 2) of macrofaunal composition of sandy beaches of Caribbean and Pacific coasts of Costa Rica. The lines show the direction of increase of the environmental data, the tidal level (I: low tide and V: high tide). PCNMs lines are predictors representing spatial eigenvectors based on geographical coordinates of sampled beaches. Proportional constrained inertia: 0.3242, Proportional unconstrained inertia: 0.6758. Only Axis 1 and 2 have p ≤ 0.001, Axis 1: Eigenvalue = 0.079, F = 22.64, Axis 2: Eigenvalue = 0.044, F = 10.79.

Fig. 6. Venn diagram for the variance partitioning (in Bold) of the effects of intertidal characteristics, geographical distance between beaches and sediments factors on the community structure. a. r2 = Adjusted r square. * = p < 0.001. 103

Estuarine, Coastal and Shelf Science 223 (2019) 94–104

J.A. Sibaja-Cordero, et al.

Acknowledgements

Dexter, D.M., 1974. Sandy beach fauna of the Pacific and Atlantic coasts of Costa Rica and Colombia. Rev. Biol. Trop. 22, 51–66. Dexter, D.M., 1976. The sandy beach fauna of Mexico. SW. Nat. 20, 479–485. Dexter, D.M., 1979. Community structure and seasonal variation in intertidal Panamanian sandy beaches. Estuar. Coast Mar. Sci. 9, 543–558. Dittmann, S., Vargas, J.A., 2001. Tropical tidal flat benthos compared between Australia and Central America. In: Reise, K. (Ed.), Ecological Comparisons of Sedimentary Shores, Ecological Studies 151. Springer, Berlin, Germany, pp. 275–293. Dray, S., Legendre, P., Peres-Neto, P.R., 2006. Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol. Model. 196, 483–493. Fauchald, K., 1977. The polychaete worms. Definitions and keys to the orders, families and genera. Nat. Hist. Mus. L.A. County, Sci. Ser. 28, 1–188. Fischer, R., 1981. Bioerosion of basalt of the Pacific coast of Costa Rica. Senckenbergiana mar 13, 1–41. Pacífico centro oriental. In: In: Fischer, W., Krupp, F., Schneider, W., Sommer, C., Carpenter, K.E., Niem, V.H. (Eds.), Guía FAO para la identificación de especies para los fines de la pesca. Pacífico Centro-Oriental, vol. I Palntas e Invertebrados. FAO, Roma. Folk, R.L., Ward, W.C., 1957. Brazos River bar: a study in the significance of grain size parameters. J. Sediment. Petrol. 27, 3–26. Fowler, L.E., 1979. Hatching success and nest predation in the green sea turtle, Chelonia mydas, at Tortuguero, Costa Rica. Ecology 60, 946–955. Fox, J., 2002. An R and S-Plus Companion to Applied Regression. Thousand Oaks, , California. Fox, J., Monette, G., 1992. Generalized collinearity diagnostics. J. Am. Stat. Assoc. 87 (417), 178–183. Frank, J.H., Ahn, K.-J., 2011. Coastal Staphylinidae (Coleoptera): a worldwide checklist, biogeography and natural history. ZooKeys 107, 1–98. Gilbert, E.R., Camargo, M.G., Sandrini-Neto, L., 2012. Rysgran: Grain Size Analysis, Textural Classifications and Distribution of Unconsolidated Sediments. R package version 2.0 [computer program]. Glynn, P.W., Dexter, D.M., Bowman, T.E., 1975. Excirolana braziliensis, a Pan-American sand beach isopod: taxonomic status, zonation and distribution. J. Zool. 175, 509–521. Hammer, Ø., Harper, D.A.T., Ryan, P.D., 2001. Past, paleontological statistics software package for education and data analysis. Palaeontol. Electron. 4, 1–9. Holme, N.A., McIntyre, A.D., 2005. Methods for the Study of Marine Benthos. Blackwells, Oxford. Jumars, P.A., Dorgan, K.M., Lindsay, S.M., 2015. Diet of worms emended: an update of polychaete feeding guilds. Ann. Rev. Mar. Sci. 7, 497–520. Keen, A.M., 1971. Sea Shells of Tropical West America. Marine Mollusks from Baja California to Peru, 2 ed. Stanford University Press, California. Krebs, C.J., 1999. Ecological Methodology, 2 ed. Addison-Welsey, California. Legendre, P., Gallagher, E.D., 2001. Ecologically meaningful transformations for ordination of species data. Oecologia 129 (2), 271–280. Lizano, O., 2006. Algunas características de las mareas en la costa Pacífica y Caribe de Centroamérica. Cienc. Tecnol. 24, 51–64. McArthur, M.A., Brooke, B.P., Przeslawski, R., Ryan, D.A., Lucieer, V.L., Nichol, S., McCallum, A.W., Mellin, C., Cresswell, I.D., Radke, L.C., 2010. On the use of abiotic surrogates to describe marine benthic biodiversity. Estuar. Coast Shelf Sci. 88, 21–32. Mourglia, V., González-Vainer, P., Defeo, O., 2015. Distributional patterns in an insect community inhabiting a sandy beach of Uruguay. Estuar. Coast Shelf Sci. 166, 65–73. O'Brien, R.M., 2007. A caution regarding rules of thumb for variance inflation factors. Qual. Quantity 41, 673–690. Ortega-Cisneros, K., de Lecea, A.M., Smit, A.J., Schoeman, D.S., 2017. Resource utilization and trophic niche width in sandy beach macro benthos from an oligotrophic coast. Estuar. Coast Shelf Sci. 184, 115–125. https://doi.org/10.1016/j.ecss.2016.11. 011. Pec, G.J., Karst, J., Taylor, D.L., Cigan, P.W., Erbilgin, N., Cooke, J.E.K., Simard, S.W., Cahill, J.F., 2017. Change in soil fungal community structure driven by a decline in ectomycorrhizal fungi following a mountain pine beetle (Dendroctonus ponderosae) outbreak. New Phytol. 213, 864–873. Powell, J., 2018. Multivariate Statistics. Western Sidney University, Hawkesbury Institute for the Environment, Sidney. Rodil, I.F., Lucena-Moya, P., Lastra, M., 2018. The importance of environmental and spatial factors in the metacommunity dynamics of exposed sandy beach benthic invertebrates. Estuar. Coasts 41, 206–217. Sibaja-Cordero, J.A., Camacho-García, Y.E., Vargas-Castillo, R., 2014. Riqueza de especies de invertebrados en playas de arena y costas rocosas del Pacífico Norte de Costa Rica. Rev. Biol. Trop. 62 (Suppl. 4), 63–84. Sibaja-Cordero, J.A., 2018. Spatial distribution of macrofauna within a sandy beach on the Caribbean coast of Costa Rica. Rev. Biol. Trop. 66 (Suppl. 1), S176–S186. Troost, A.I., Rupert, S.D., Cyrus, A.Z., Paladino, F.V., Dattilo, B.F., Peters, W.S., 2012. What can we learn from confusing Olivella columellaris and O. semistriata (Olivellidae, Gastropoda), two key species in panamic sandy beach ecosystems? Biota Neotropica. https://doi.org/10.1590/S1676-06032012000200011. Villacampa, Y., López, I., Aragonés, L., García, C., López, M., Palazón, A., 2017. Water quality of the beach in an urban and not urban environment. Int. J. Sustain. Dev. Plan. 12 (4), 713–723.

This study is part of project 808-B4-117 of CIMAR, partially funded by Vicerrectoría de Investigación, Universidad de Costa Rica and and the TWAS/CONICIT award of The World Academy of Sciences, Italy and Consejo Nacional para Investigaciones Científicas y Tecnológicas, Costa Rica.. We are grateful to the students of the School of Biology, Universidad de Costa Rica S. Pastor, N. Goebel, B. Antillón, R. Barboza, J. López, O. Segura, J.M. Valverde, C. Salas, R. Cambronero, S. Orozco and V. Castellví, who helped at processing samples in the laboratory of CIMAR. Paul Hanson of the School of Biology, Universidad de Costa Rica helped with the insect identification. David Butvill helped with the revision of the English. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.ecss.2019.04.032. References Alongi, D.M., 1989. Ecology of tropical soft-bottom benthos: a review with emphasis on emerging concepts. Rev. Biol. Trop. 37, 85–100. Alongi, D.M., 1990. The ecology of tropical soft-bottom benthic ecosystems. Oceanogr. Mar. Biol. Annu. Rev. 28, 381–496. Barboza, F.R., Defeo, O., 2015. Global diversity patterns in sandy beach macrofauna: a biogeographic analysis. Sci. Rep. 5, 14515. Borcard, D., Legendre, P., Avois-Jacquet, C., Tuomisto, H., 2004. Dissecting the spatial structure of ecological data at multiple scales. Ecology 85 (7), 1826–1832. Brusca, R.C., Iverson, E.W., 1985. A guide to the marine isopod Crustacea of Pacific Costa Rica. Rev. Biol. Trop. 33 (Suppl. 1), 1–77. Celentano, E., Lercari, D., Maneyro, P., Rodríguez, P., Gianelli, I., Ortega, L., Orlando, L., Defeo, O., 2019. The forgotten dimension in sandy beach ecology: vertical distribution of the macrofauna and its environment. Estuar. Coast Shelf Sci. 217, 165–172. https://doi.org/10.1016/j.ecss.2018.11.008. Checon, H.H., Corte, G.N., Shah Esmaeili, Y.M.L., Amaral, A.C.Z., 2018. Nestedness patterns and the role of morphodynamics and spatial distance on sandy beach fauna: ecological hypotheses and conservation strategies. Sci. Rep. https://doi.org/10. 1038/s41598-018-22158-3. Coan, E.V., Valentich-Scott, P., 2012. Bivalve seashells of tropical west America. Marine Bivalve Mollusks from Baja California to Northern Perú. Santa Barbara Museum of Natural History Monographs 6, Studies in Biodiversity 4, Santa Barbara, California. Colombini, I., Bouslama, M.F., Elgtari, M., Fallaci, M., Scapini, F., Chelazzi, L., 2005. Study of the community structure of terrestrial arthropods of a Mediterranean sandy beach ecosystem of Morocco. In: Bayed, A., Scapini, F. (Eds.), Ecosystèmes Côtiers Sensibles de La Méditerrané: Cas Du Littoral de Smir (Série Générale, 4). Travaux del' Institut Scientifique, Rabat, pp. 43–54. Corrales-Ugalde, M., Sibaja-Cordero, J.A., 2015. Macrofauna bentónica de las playas de arena del Área de Conservación Osa, Puntarenas, Pacífico Sur de Costa Rica. Rev. Biol. Trop. 63 (Suppl. 1), 273–285. Cortés, J., 2014. Compilación y análisis de las investigaciones científicas sobre temas marinos y atmosféricos en el Pacífico Norte de Costa Rica. Rev. Biol. Trop. 62 (Suppl. 4), 151–184. Cortés, J., Fonseca, A.C., Nivia-Ruiz, J., Nielsen-Muñoz, V., Samper-Villarreal, J., Salas, E., Martínez, S., Zamora-Trejos, P., 2010. Monitoring coral reefs, seagrasses and mangroves in Costa Rica (CARICOMP). Rev. Biol. Trop. 58 (Suppl. 3), 1–22. Crawley, M.J., 2007. The R Book. John Wiley & Sons, Chichester. de León-González, J.A., Bastida-Zavala, J.R., Carrera-Parra, L.F., García-Garza, M.E., Peña-Rivera, A., Salazar-Vallejo, S.I., Solís-Weiss, V. (Eds.), 2009. Poliquetos (Annelida: Polychaeta) de México y America Tropical. Universidad Autónoma de Nuevo León, Monterrey. Dean, W.E., 1974. Determination of carbonate and organic matter in calcareous sediments and sedimentary rocks by loss on ignition; comparison with other methods. J. Sediment. Res. 44, 242–248. Defeo, O., McLachlan, A., 2013. Global patterns in sandy beach macrofauna: species richness, abundance, biomass and body size. Geomorphology 199, 106–114. Defeo, O., Barboza, C.A.M., Barboza, F.R., Aeberhard, W.H., Cabrini, T.M.B., Cardoso, R.S., Ortega, L., Skinner, V., Worm, B., 2017. Aggregate patterns of macrofaunal diversity: an interocean comparison. Glob. Ecol. Biogeogr. 26, 823–834. Defeo, O., McLachlan, A., Schoeman, D.S., Schlacher, T.A., Dugan, J., Jones, A., Lastra, M., Scapini, F., 2009. Threats to sandy beach ecosystems: a review. Estuar. Coast Shelf Sci. 81, 1–12. Dexter, D.M., 1972. Comparison of the community structures in a Pacific and an Atlantic Panamanian sandy beach. Bull. Mar. Sci. 22, 449–462.

104