Sandy Beach Macrofauna of Yucatán State (Mexico) and Oil Industry Development in the Gulf of Mexico: First Approach for Detecting Environmental Impacts

Guerra-Castro, Edlin; Hidalgo, Gema; Castillo-Cupul, Raúl E.; Muciño-Reyes, María; Noreña-Barroso, Elsa; Quiroz-Deaquino, Jaime; Mascaro, Maite; Simoes, Nuno

doi:10.3389/fmars.2020.589656

ORIGINAL RESEARCH article

Front. Mar. Sci., 29 October 2020

Sec. Marine Ecosystem Ecology

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.589656

This article is part of the Research Topic Marine Biodiversity Observation Network (MBON) View all 16 articles

Sandy Beach Macrofauna of Yucatán State (Mexico) and Oil Industry Development in the Gulf of Mexico: First Approach for Detecting Environmental Impacts

$\r\nEdlin Guerra-Castro,*&#x;$ Edlin Guerra-Castro^1,2*†

Gema Hidalgo^2,3†

Raúl E. Castillo-Cupul^4†

María Muciño-Reyes^5,6†

Elsa Noreña-Barroso^2,7†

Jaime Quiroz-Deaquino⁵

Maite Mascaro^2,5†

Nuno Simoes^2,5,8*†

¹Escuela Nacional de Estudios Superiores Unidad Mérida, Universidad Nacional Autónoma de México, Mérida, Mexico
²Laboratorio Nacional de Resiliencia Costera, Laboratorios Nacionales, CONACYT, Mexico City, Mexico
³Centro del Cambio Global y la Sustentabilidad A.C., Villahermosa, Mexico
⁴Kalanbio A.C., Mérida, Mexico
⁵Unidad Multidisciplinaria de Docencia e Investigación Sisal (UMDI-SISAL), Facultad de Ciencias, Universidad Nacional Autónoma de México, Sisal, Mexico
⁶Posgrado en Ciencias Biológicas, Universidad Nacional Autónoma de México, Ciudad de México, Mexico
⁷Unidad de Química Sisal, Facultad de Química, Universidad Nacional Autónoma de México, Sisal, Mexico
⁸International Chair for Coastal and Marine Studies in Mexico, Harte Research Institute for Gulf of Mexico Studies, Texas A&M University-Corpus Christi, Corpus Christi, TX, United States

The biodiversity of the coastal ecosystems in the Gulf of Mexico is threatened by anthropogenic activities of various kinds. The predominant portion of the land-sea margin in the State of Yucatán consists of exposed sandy beaches. This ecosystem is threatened by several activities that vary in spatial scales, at a local: cargo/fishing ports, touristic facilities, maritime traffic, and domestic pollution; and at a larger scale: the forthcoming development of the oil industry. In the absence of information about the biodiversity of the beaches of Yucatán, we implement the Marine Biodiversity Observation Network (MBON) Pole to Pole sampling protocol to (1) Quantify the spatial patterns of diversity of macrofauna along the beaches; (2) quantify current levels of pollution by hydrocarbons (aromatic and aliphatic); (3) estimate sampling effort for future environmental impact assessments. During November 2018, six localities along the coastline of Yucatán State were sampled following a spatial hierarchical design that included three sites at each locality, and 9–18 core samples in the intertidal strata of each site. As a result, 31 species of invertebrates were registered. The patterns of distribution and abundance of species showed that there was a base community structure along the entire coast dominated by four species. In general, the density of species was relatively low (2–4 species/0.01 m³) and the density of individuals was high (20–200 ind./0.01 m³). The beta diversity was higher between localities with good environmental health, but the estimated alpha diversities did not show a pattern regarding the health of the coast. Lastly, the overall number of species reported suggest that gamma diversity of the macrofauna in the beaches of Yucatán is within the highest known worldwide. Levels of hydrocarbons detected in this study are exceptionally low compared to those reported for other coastal areas of the Gulf of Mexico and are several orders of magnitude lower than those considered as lower-threshold values for marine sediments. The biological and chemical patterns reported here indicated that it is an appropriate moment to start long-term monitoring. The sampling design we suggest is based on statistical precision and internationally recognized protocols for assessment of marine diversity on sandy beaches.

Introduction

Biodiversity of the coastal ecosystems in the Gulf of Mexico is threatened by anthropic activities of various kinds (Peterson et al., 1996; Andrade, 2010; Joye, 2016). Especially on the Southern and North-western coast, the environmental impacts of oil exploitation and marine traffic stand out. In 2013, the Gulf of Mexico contributed with 54% of oil and 47% of the US production, whilst 75% of Mexico oil production comes from this area (National Ocean Service NOAA, 2011; Yoskowitz et al., 2013). In the same year, there were 3,095 oil & gas related offshore platforms in the northern GMx coast (United States), and 183 in the southern GMx coast (Mexico), with many more expected to be built due to the exploration of deeper fields within the US jurisdiction, and new areas in Mexico as a result of the recent energy reform in Mexico (SENER, 2015). Therefore, there is a historic track with many oil spill accidents, including two very large ones with enormous environmental consequences: the Ixtoc I by PEMEX (Bay of Campeche, 1979; see Sun et al., 2015), and the Deepwater Horizon (Delta of the Mississippi, 2010) by BP (Jernelov and Linden, 1981; Jernelov, 2010). Currently, no oil is extracted from the Yucatán platform, but that condition will change as foreseen in the 5-year Mexican plan for the exploration and extraction of hydrocarbons 2015–2019 (SENER, 2015). Of the 12 Oil Mexican Provinces recognized by the Mexican Hydrocarbons Commission, two are prone to generate environmental impacts on the coasts of the State of Yucatán. The first is the shallow water exploration area AS1012, with a coverage of 461 km², located on the Yucatán platform, and whose reserves are estimated at 13.5 MMbpce (SENER, 2015). The second is the Deep Gulf of Mexico, which is the province with the greatest potential for exploitation of conventional resources, estimated at 27.8 MMbpce (SENER, 2015). Besides, it was recently announced the construction of a fuel terminal in the Puerto de Altura in Progreso (Figure 1), with a capacity to store 70 million liters of fuel. While all these activities will contribute to the economic development of the nation, they will also greatly increase the potential of environmental impacts of different magnitudes on regional marine biodiversity (Andrade, 2010; Pech Pool et al., 2010).

FIGURE 1

Figure 1. Localities studied on the coast of Yucatán State, México. Three portions of the coast are colored according to environmental health (sensu LANRESC, 2017): orange (regular), red (poor), green (good). For each locality, three sites distance apart by 3–9 km were sampled. Also, for each locality, a picture of the beach is shown.

The coast of the Yucatán is characterized by its calcareous origin, without surface drainage, and consists mostly of dissipative beaches of low slope (Ramos, 1975). The coastal ecosystems include sandy beaches, mangrove forests, seagrass meadows, and reef archipelagos that increase the physiographic complexity and importance of the natural capital of the region. There has been an important scientific descriptive effort of biodiversity and environmental characteristics of mangrove forests and coastal lagoons (Herrera-Silveira et al., 1998; Tapia González et al., 2008), as well as reef archipelagos of this region (González-Muñoz et al., 2013; Ortigosa et al., 2015; Ugalde et al., 2015; Mendoza-Becerril et al., 2018; Palomino-Alvarez et al., 2019; Paz-Ríos et al., 2019, 2020; Robertson et al., 2019); however, the biodiversity of the beaches has been little studied. Sandy beaches represent around 286 km (86%) of coastline in the State of Yucatán. These beaches are critical habitats for shorebirds (Arturo Lopez et al., 1989), turtles (Cuevas et al., 2010), as well as dune vegetation along the entire coast (Espejel, 1984; Islebe et al., 2015). They also support diverse touristic activities, especially important for the local economy (Meyer-Arendt, 2001; Cuevas Jiménez et al., 2016). Due to its extension, it is the one that presents the greatest probability of being negatively impacted by the potential accidents of the oil & gas industry in the Yucatán platform.

In general, impacts generated by the oil & gas industry in coastal ecosystems imply loss and long-term modification of local biodiversity (Teal and Howarth, 1984; Kingston, 2002). The effects of oil spills on the biodiversity of sandy beaches have strongly depended on the magnitude of the spill, on the granulometric properties of the beaches and the persistence of oil in the sediment (Bejarano and Michel, 2016). Estimating these effects requires a good understanding of the patterns of spatial and temporal variation that naturally occur on sandy beaches. Likewise, to identify the level of resilience of the ecosystem (i.e., biological recovery), it will be essential to know the system well before the disturbance, and this implies estimating the spatiotemporal variation of its ecological components at different spatial and temporal scales (Underwood, 1991; Cruz-Motta et al., 2007; Bejarano and Michel, 2016). The information necessary to detect environmental impacts of this nature arises from the so-called “Baseline Studies,” framed within the “Environmental Impact Assessment approach” (Mareddy, 2017). However, the vast majority of these studies have transcendental failures in the sampling designs, some of them are: (a) the measured impact indicator variables are inadequate (e.g., Shanon index of diversity), (b) the sampling start shortly before the project execution process (e.g., few weeks before the start of the operations), (c) there are not reference localities and, (d) the sample size is not representative of the natural variability (e.g., usually three-five samples per site), resulting in statistical tests of low power (Underwood and Chapman, 2003). Consequently, there is rarely enough solid information to unequivocally indicate the existence of an environmental impact (Bulleri et al., 2007; Cruz-Motta et al., 2007), this is particularly remarkable in post-impact environmental assessments on polluted sandy beaches (Schlacher et al., 2008; Bejarano and Michel, 2016).

Sandy beaches provide several ecosystem services to society, but they are threatened by various human activities (Schlacher et al., 2008; Defeo et al., 2009; Rodríguez-Revelo et al., 2018; Martínez et al., 2020). In general, some of the recognized anthropogenic threats with the greatest risk on beaches are surface physiographic removal, microbial biohazards, introduced technological hazards, chronic chemical hazards, and chronic geopolitical hazards (Fanini et al., 2020). The risk associated with oil spills is considered lower than those mentioned (Fanini et al., 2020), although they have the potential to generate mass mortalities and drastically reduce the quality of habitats for several years (Bejarano and Michel, 2016). Unfortunately, this habitat has been ignored in long term ecological surveys. There is a global need to measure beaches’ resilience to short (pulse) and long (press) term perturbations and increased demand for information from less-represented areas, including beach morphodynamics, pollution, ecological and socio-economic indicators, to undertake comprehensive and long-term impact assessments (Fanini et al., 2020; Thom, 2020).

Although there are significant efforts to compare spatiotemporal patterns of diversity in these ecosystems on a global scale (Defeo and McLachlan, 2005, 2013; McLachlan and Dorvlo, 2005; Rodil et al., 2014; Barboza and Defeo, 2015), a lack of standardized methodologies for this very dynamic environment, with huge variation in geological history, tide range, sediment texture, slope, and exposition to waves, make comparisons of spatial patterns a big challenge, compromising the forecasting of future ecological scenarios and management activities. To overcome this limitation, some workshops have been held to develop sampling protocols that allow the comparison of results in regional contexts (Schlacher et al., 2008; Canonico et al., 2019). To attend this, the Marine Biodiversity Observation Network - Pole to Pole of the Americas (MBON Pole to Pole) promoted the design of a standardized sampling protocol for sandy beaches that allows obtaining comparable information on all continental coasts (MBON Pole to Pole, Pole, 2019). This protocol was recognized as an OceanBestPractice initiative of UNESCO/IODE¹.

Considering the imminent development of oil and gas exploration and production activities on the Yucatán platform, as well as the background of environmental impact in other nearby regions of the Gulf of Mexico, it is predictable the potential generation of negative impacts on the biodiversity of the coastal ecosystems of the State of Yucatán, mainly its beaches. In this sense, a baseline study with biological indicators at different spatial scales, which allows detecting non-natural changes in the biological component, as well as the potential chemical mechanisms that may cause any loss or alteration of biodiversity, are required. To help address this need, this study presents an extensive pilot data, with biological and chemical variables, as a first baseline data set of the region, as well as a quantitative description of the patterns of diversity of macrofauna and pollutants in the Yucatán Peninsula, using the MBON Pole to Pole sampling protocol. Finally, we make recommendations regarding the sampling effort necessary for future environmental assessment.

Materials and Methods

Study Area

Sampling was carried out in November 2018 on six localities (distance apart by 50 km or more) of the windward side of the Yucatán peninsula; from west to east: Celestún, Sisal, Progreso, Telchac, Dzilam and El Cuyo (Table 1 and Figure 1). The entire littoral is characterized by being microtidal, sea breeze dominated with a dissipative profile, low slopes, and medium to coarse sands (Appendini et al., 2012; Medellín and Torres-Freyermuth, 2019). In general, the region is karstic, with no superficial rivers or estuaries, but with well-developed mangrove forests surrounding coastal lagoons that are strongly influenced by groundwater discharges (Herrera-Silveira et al., 1998; Tapia González et al., 2008). These coastal lagoons are connected with the sea in several zones along the littoral, being the source of hipohaline waters with a high concentration of dissolved inorganic nutrients, especially in the rainy season. Trade winds tend to dominate in spring and summer, but in autumn and winter, the pattern changes with recurrent cold fronts from the north, with increased macrophyte wrack deposition (Enriquez et al., 2010; Medellín and Torres-Freyermuth, 2019). The current environmental health of the coast was assessed in 2017 by an interdisciplinary panel (LANRESC, 2017), resulting in three heterogeneous coastal segments. The portion of the coast between Celestún and Sisal is recognized as an area with regular environmental health; then, beyond Sisal until Dzilam, including Progreso and Telchac, it is a region classified with areas with poor environmental health; finally, Dzilam and El Cuyo stand out for presenting good environmental health (Table 1 and Figure 1).

TABLE 1

Table 1. Names of localities sampled along the coast of Yucatán.

Sampling Protocol

The sampling design included three spatial scales: tens of kilometers (localities), few kilometers (sites), few meters (cores). At each locality, three sites (separated by 3–9 kilometers) were sampled. None of the sites was less than 2 km from the mouth of any coastal lagoon or estuary. We applied the first version of the MBON P2P sampling methodology for sandy beaches, inspired in recommendations by Schlacher et al. (2008). At each site (a.k.a. sampling station), three transects (separated by 10 m) were displayed perpendicularly to the coastline. Along each transect, a sample of sand was collected each meter, starting from the waterline to the drift line (measured in all cases with a tape measure), this protocol ensures collecting samples from the infra, meso, and supralittoral. The fauna was collected with a core of 22 cm in diameter, which was inserted between 20 and 25 cm inside the sediment, representing 0.038 m² of area or 0.01 m³ of volume. The tidal strata (infra, meso, supra), as well as the depth of the core, were annotated for each sample core. The samples were transferred fresh to the laboratory for immediate processing. In addition to biological samples, four sediment samples were taken at one randomly chosen site of each locality to estimate the concentration of aliphatic and aromatic hydrocarbons, as well as for granulometric analysis. These samples were taken at the meso littoral zone using the same core used for biological sampling.

Laboratory Protocols

Biological Samples

Samples were sieved with a 0.5 mm mesh opening, and the retained organisms were photographed alive and subsequently fixed in a 4% solution of neutralized formaldehyde with sodium tetraborate. Organisms were identified to the lowest possible taxonomic level and counted. The specimens were deposited in the biological collection of the UMDI-Sisal, UNAM. Also, a database was built using the standard Darwin Core (Darwin Core Task Group, 2009)². This database is available from the Global Biodiversity Information Facility³ and visible on the Ocean Biodiversity Observation Network⁴. Alternatively, all the data can be accessed at Zenodo (DOI: 10.5281/zenodo.3771828) (Guerra-Castro et al., 2019).

Sand Samples

Sediment samples collected for hydrocarbon determination were kept frozen at −20°C and then freeze-dried and passed through a 500 μm sieve. Concentrations of aliphatic and aromatic hydrocarbons in sediment (15 of the 16 PAHs considered as priorities by the US Environmental Protection Agency, EPA) were determined based on modifications of EPA methods 3550C, 3535 and 8270D (US EPA, 2007a,b,c). For each sample, 6 g of lyophilized sediment were extracted twice by ultrasound-assisted extraction (USE) with 12 mL hexane:acetone (1:1, v/v) using an ultrasonic processor (Cole Palmer CPX500) at 60% amplitude over 2 min. After each extraction, the organic phase was separated by centrifugation (5,000 rpm for 10 min). Extracts were mixed, treated with activated copper to remove sulfur, and concentrated using a rotary-evaporator. Hydrocarbon fractions were obtained from extracts by solid-phase extraction (SPE) using C-18 500 mg/6mL cartridges (Supelclean ENVI-18, 57064, Supelco). Cartridges were conditioned with 10 mL of hexane, sample extracts were passed by gravity flow and then eluted with 10 mL of hexane to obtain the aliphatic fraction, and 5 mL of hexane:dichloromethane (7:3, v/v) followed by 5 mL of dichloromethane to recover the PAH fraction. Fractions were evaporated using a gentle nitrogen flow and individual hydrocarbons were determined by gas chromatography/mass spectrometry (GC-MS), using a gas chromatograph coupled to a mass selective detector operated in electron impact (EI) ionization mode and equipped with an automatic liquid sampler (Agilent Technologies 7890B Series GC; 5977B MSD and 7693A Autoinjector, respectively). The injection was carried out in split-less mode (1 min) at 280°C. Chromatographic separation was performed using a J&W HP-5MS capillary column (30 m × 0.25 mm I.D. and 0.25 μm of film thickness). Carrier gas was He (ultra-pure grade) with a flow rate of 0.8 mL/min; oven temperature was initially set at 60°C, then increased 6°C/min to 290°C (hold time 11.67 min). Mass spectra (m/z 50–550) were recorded at a rate of five scans per second at 70 eV. Mass spectrometric analysis for quantitative determination was performed by selected ion monitoring (SIM Mode) of two characteristic fragment ions for each analyte (Supplementary Table 1). Analytical quality control included procedural blanks, calibration curves, and internal standards. The aliphatic hydrocarbon detection limit was 0.5 ng/g and PAH detection limits ranged from 0.09 to 0.79 ng/g. Particle size analyzes were done with conventional sieving methods (Keith, 1996).

Analyses of Data

Patterns of Species Diversity

We evaluated patterns of spatial variation for four features of intertidal macrofauna diversity: (1) structure of the assemblages, (2) patterns of density of species and abundance of organisms (3) beta diversity, and (4) alfa and gamma diversity. For the structure of the assemblages, the counts of each species were arranged in an N × P matrix, with N being the total number of samples and P being the number of species. This matrix was transformed into the natural logarithms (plus 1) to downweigh the effect of highly abundant species. Then, the Bray-Curtis coefficient of dissimilarity was estimated between each pair of samples, generating a matrix of dissimilarities that was used for statistical analyses and ordination. The total variation of this matrix was decomposed with a multifactorial linear model using Permutational Multivariate Analysis of Variance (Anderson, 2017). In the model, the main sources of variation were: Localities (fixed factor with six levels: Celestún, Sisal, Progreso, Telchac, Dzilam, and El Cuyo), Sites (random factor, nested in localities, three levels: site 1, site 2, site 3) and Strata (fixed factor, three levels: infra, meso, supra). The null hypotheses were generated using 9,999 permutations of residuals under the reduced model (Anderson and Ter Braak, 2003). The relative importance of each spatial scale (tens of km, few km, few m) were identified as the relativized square root of the pseudo-component of variation of Localities, Sites, and residuals in the model. Then, to visualize the pattern of similarity among localities, centroids for each locality-site-strata were estimated and projected with a non-metric MDS. Patterns of distribution and abundance of species along the coast were identified with an ordered shade plot using constrained seriation of most similar species in its spatial distribution regarding the locality (Clarke et al., 2014b). Variability in the density of species and abundance (i.e., number of species per sampling core, and the number of individuals per sample core, respectively), were analyzed using the same linear model for the structure of the assemblage but with univariate analysis of variance based on permutations. Both variables were analyzed without any transformation and using 9,999 permutations of residuals under the reduced model. Patterns of spatial trends in the density of species and abundance were summarized and represented using mean (±SD) plots. All these analyses were done using the software PRIMER v7 & PERMANOVA (Clarke et al., 2014a).

For beta diversity, we explore the relationship between the Jaccard dissimilarities in species composition along the coast of Yucatán (approach T3 sensu Anderson et al., 2010) and the partitioning proposed by Baselga (2010). To identify the proportion of species replacement and species loss along the coast, the Jaccard dissimilarities among sites were estimated and partitioned in their components of turnover/nestedness and correlated with the geographic distances among sites using the rank-based Spearman coefficient of correlation. The null hypothesis of no spatial autocorrelation was tested using 9,999 permutations. These analyses were done with the statistical software R (R Core Team, 2013) using the packages betapart (Baselga and Orme, 2012) and vegan (Oksanen et al., 2015).

Alfa diversity was associated with the potential number of species per locality, while gamma diversity as the potential number of species along the coast of Yucatán. In both cases, we use species accumulation curves with interpolation-extrapolation of Hill Numbers and the incidence-based non-parametric estimator Chao2 (Chao and Jost, 2012). Considering that the number of sampling was not the same in all localities, we also estimate the sample completeness for each one, this measure allows us to infer the representativeness of the sampling effort at each locality. The interpolation of completeness was used to compare the local richness for a common completeness value of 0.9 (Chao et al., 2009; Chao and Jost, 2012). The analyses were done with the statistical software R (R Core Team, 2013) using the package iNEXT (Hsieh et al., 2016).

Patterns of Hydrocarbons Pollution and Environmental Variation

Individual concentrations that were less than the limits of detection were assigned a value of zero for statistical estimations. Concentrations of PAHs were classified and summed according to the number of rings of each molecule (i.e., 2–3 rings, 4 rings, 5–6 rings), and then totalized (ΣPAHs). Similarly, aliphatic hydrocarbons were classified and summed according to the molecular structure (odd and even carbon number) and summed as total (Σn-alkanes). To identify the potential source of aliphatic compounds, the Carbon Preference Index (CPI) was estimated as:

\begin{matrix} \frac{\sum odd carbon number n - alkanes}{\sum even carbon number n - alkanes} \end{matrix}

where values greater than 1 indicate natural sources, while values lesser than 1 indicate anthropogenic sources (Marzi et al., 1993). Then, patterns of spatial heterogeneity in PAHs and n-alkanes were visualized using standard statistical representations of totals and the respective fraction of classified molecules. Besides, the null hypothesis of equal concentrations of these molecules between localities was tested using multivariate analyses. For this, data was transformed to a common scale (centered in mean zero with variance unit, aka normalization), then, a matrix of Euclidean distances was estimated, and total variation partitioned using a one-way permutational analysis of variances. For this, 9999 permutations of raw data were used to generate the null distribution of pseudo-F values.

Sampling Effort for Future Assessments

The criterion for redefining sampling effort was the optimization of precision, that is, estimating the number of samples to improve the precision in the estimates of variability in assemblage structure and concentration of hydrocarbons at a reasonable cost, in terms of the number of samples at each locality. For this, we use data of this study as pilot data. For changes in the structure of the assemblages, multivariate standard error (MultSE) (Anderson and Santana-Garcon, 2015) were analyzed for several sampling efforts using simulations of data and the R package SSP (Guerra-Castro et al., 2020). For each stratum at each locality, 20 virtual sites were simulated, each one with N = 100 potential sample cores. Resampling were done for each combination of n = 2 to n = 20 and sites m = 2 to m = 20. For each combination of n and m, MultSE was estimated; this estimation was repeated 10 times. All these processes were repeated for 10 simulated data sets. Then, a metanalysis about the behavior MultSE for each sampling effort was done projecting the MultSe – sampling effort relationship. The optimal sampling effort was defined as the range in which an additional sampling unit improves the worst precision by 10%, but not beyond 2.5%. Beyond this point, it was considered unnecessary sampling efforts. For pollutants, we used the equation $n = {[s / (p \bar{x})]}^{2}$ , being n the sampling effort, s the standard deviation, $\bar{x}$ de arithmetic mean of the sample, and p a standardized precision (i.e., 0.2) (Andrew and Mapstone, 1987). The variables used for these estimations were ΣPAHs and Σn-alkanes, and the estimations were done for each locality. The codes for these analyses are available as Supplementary Material.

Results

Patterns of Species Diversity

A total of 31 taxa/species were registered in 225 core samples, of which 15 were arthropods (13 crustaceans), 10 annelids (mainly polychaetes), 4 mollusks, and 2 bryozoans. The list of species (including main biological traits) is available in Supplementary Table 2, and the spatial distribution of each species can be visualized at OBIS (Guerra-Castro, 2019) and accessed at Zenodo (Guerra-Castro et al., 2019). Significant statistical differences in the structure of the assemblage were detected in all the spatial scales (Table 2A, p-values <0.05). The potential effect of the zone (i.e., littoral strata) varied among sites (Table 2A, interaction term Strata x Site (Locality) with p-value <0.05). The largest source of variation was the residuals: cores of the same littoral strata in the same site. This source of variation was about 34% in Bray Curtis dissimilarity and about 44% of the total variation. The second source of variation was the difference among localities (square root of components of variation 22% of Bray-Curtis dissimilarity), followed by variability among sites (18%) and difference among littoral strata (13%). El Cuyo and Dzilam were the most different localities (Figure 2), the rest of the localities were more similar between them but different from El Cuyo and Dzilam (Figure 2). Despite these differences, it can be noted that all sites present a community structure with the same four dominant species: the polychaete Polyophthalmus pictus, the oligochaete Tubificoides diazi, and the isopods Excirolana mayana and E. braziliensis (Figure 3). In general, differences among localities were explained by a pool of 4–7 species confined to each locality (Figure 3). Regarding beta diversity, the total multi-site Jaccard dissimilarity was 0.91, indicating that only 9% of species were shared across all sites. In general, the turnover component was considerably higher (0.87) than the nestedness component (0.04). This spatial pattern of dissimilarities was weakly correlated with spatial distances (ρ = 0.26, p < 0.001).

TABLE 2

Table 2. Permutational Analyses of Variances to evaluate spatial patterns of variability in the structure of intertidal macrofauna in sandy beaches of Yucatán, according to a multifactorial linear model that includes locality (L), site (S) and littoral strata (St).

FIGURE 2

Figure 2. Non-metric Multidimensional Scaling of centroids for the interaction term Strata × Site (Locality), of macrobenthic organisms in sandy beaches of Yucatán. The abundances of species were log-transformed. The coefficient of dissimilarity used was Bray-Curtis. Numbers over each symbol represent the site ID.

FIGURE 3

Figure 3. Shade plot of the averaged samples (columns) according to the locality, site (1, 2, 3), and littoral strata (triangle for supra, square for meso, circle for infra), ordered from West to East of the Yucatán coast. The depth of gray shading is linearly proportional to a logarithmic transformation (+1) of the counts of species (e.g., e^5.6 – 1 = 269 ind., e^3.5 – 1 = 32 ind.). The order of species (rows) corresponds to the similarity of the distribution patterns of the species, based on a hierarchically agglomerative clustering using Whittaker’s index of Association over standardized abundances.

The density of species, as well as the density of individuals, were statistically significant at the scale of locality (Tables 2B,C, factor Locality with p-values <0.05). These metrics were higher in Dzilam (Figure 4), specifically, the average richness per core was around 4 species, while the density was very variable but always over 100 individuals, with a mean around 200 individuals. In the other localities, the density of species averaged between 2 and 3, and the number of individuals averaged between 25 and 50 individuals. The effect of strata was significant for the density of species, but such effect was variable among sites of the same locality (Table 2B, interaction term Strata × Site (Locality) with p-value <0.05). No statistical variability in the density of individuals was detected among sites, neither difference between littoral strata.

FIGURE 4

Figure 4. Overall means and standard errors for the density of species (A) and density of organisms (B) per core across the sandy beaches of Yucatán. The standard error was calculated as the standard deviation of the Site averages of each locality. Localities are presented from West to East.

The localities with higher richness registered were Dzilam and El Cuyo, with 15 species each. Paradoxically, Dzilam was the locality with the smallest number of samples (26 cores, due to the narrow littoral zone), while El Cuyo was the locality with a larger sampling effort (44 cores). On the other hand, the localities with the lower richness of species were Celestún and Progreso, both with 9 species (Figure 5A). Between these extremes of diversity, Telchac and Sisal account for 12 and 13 species, respectively (Figure 5A). In all cases, the sample completeness was greater than 0.9 (Figure 5B), which supposes a deficit in the detection of species lesser than 10% in all localities. After the evaluation of the interpolation of richness for completeness at 0.9, it was detected that Dzilam, El Cuyo, Sisal, and Telchac are the localities with higher richness, while Celestún and Progreso were the sites with lower species richness (Figure 5C). The gamma diversity, extrapolated using the Chao2 estimator, could be greater than 60 species. The actual sample completeness is 0.885, however, from Figure 5D it seems that the richness of species of the Yucatán beaches could be much higher than that reported in this study.

FIGURE 5

Figure 5. Interpolation and extrapolation of richness values and sample completeness analyses of sandy beaches of Yucatán. (A) Species accumulation curve based on the incidence of species for each locality. (B) Sample-coverage accumulation curve based on incidence for each locality. (C) Sample completeness curves linking curves in panels (A,B). (D) Species accumulation curve for the coast of Yucatán based on the incidence of species in 18 sampled sites. The solid lines represent interpolation, while dashed lines extrapolation. Solid circle represents the observed value.

Patterns of Hydrocarbons Pollution and Granulometry

Average and standard deviation of the individual concentration of 16 PAH’s and 28 aliphatic hydrocarbons are listed in Supplementary Table 3. In general, the highest values of ΣPAHs were measured in Progreso, followed by Sisal and Telchac (Figure 6A). Molecules with 2–3 rings tend to predominate in all localities except in Progreso, where the concentrations of 4 rings and 5–6 rings molecules were as high as 2–3 ring PAH content. PAHs with 4 rings and 5–6 rings were 2–3 times higher in Progreso than in other localities (Figure 6A). Similarly, the highest values of Σn-alkanes were measured in Sisal, Progreso, and Telchac, doubling values concerning other localities (Figure 6B). In all cases, the CPI was around to unity (Figure 6B). Besides the apparent spatial trends in PAH’s and aliphatic hydrocarbons along the coast of Yucatán, especially nearby Progreso (Figures 6A,B), there was no statistical evidence to reject the null hypothesis of equal concentration of pollutants across these localities (PERMANOVA, Pseudo-F = 1.38, p-value >0.05). Although it is not possible to perform an analysis of power on this test to assess the probability of type II error (failure to reject the null hypothesis), the low number of samples per location (n = 4) and the high variability within locations may have influenced the significance of the test.

FIGURE 6

Figure 6. Analyses of sediments on the beaches of the state of Yucatán. (A) Average and standard deviation of the Total of PAHs and their respective fractions according to the number of rings. (B) Average and standard deviation of the Total of Aliphatic Hydrocarbons, and their fractions according to the number of carbon atoms. The CPI (mean ± sd) for each locality is showed over each column. (C) Average of the relative fractions of the grain size of the sediments.

Regarding granulometry, a statistically significant difference was found between localities (PERMANOVA, Pseudo-F = 7.23, p-value <0.05); specifically, El Cuyo compared to the other localities (Pair-wise t-test, p < 0.05). Except for El Cuyo, the entire coast was characterized by sediments with a greater proportion of medium to coarse and very coarse sands (>80%, Figure 6C). El Cuyo was characterized by having more abundant fine and very fine sands fractions (>70%, Figure 6C) than in the rest of the localities (<20%, Figure 6C).

Estimation of Sampling Effort for Future Assessments

The optimal sampling effort to characterize the macrofauna ranged between 7 and 10 samples per strata (Figure 7, dark gray area in upper panels). This sampling effort would allow obtaining a MultSE between 0.9 and 0.175, corresponding the 45 and 55% of the MultSE obtained with a sampling effort of two samples and two sites, respectively. The analysis highlights that the estimated range of effort varies according to stratum and localities. Consistently, Dzilam is the site that requires the lower sampling effort in all strata; on the other hand, El Cuyo requires the greatest sampling effort for the supra and infralittoral strata, and Telchac the locality that requires the greatest sampling effort in the mesolittoral strata (Figure 7, dark gray area in upper panels). In all cases, the sampling range should be between 7 and 10 samples per strata at each site. Regarding the number of sites, the optimal improvement number was estimated between 4 and 11 sites per locality. This effort would reduce the MultSE from 0.03 to 0.009–0.02 (Figure 7, dark gray area in lower panels). Based on these estimates, we recommend increasing the number of samples for strata to 8, and the number of sites to 4. This recommendation will be analyzed in section “Discussion.”

FIGURE 7

Figure 7. Behaviors of MultSE for each stratum in each of the localities. 10 simulations were performed, in each one 20 virtual sites were generated for each stratum, with 100 potential sampling units. In each simulation, n = 2 to n = 20 and sites m = 2 to m = 20 were resampled, each combination of n and m was repeated 10 times. For each, the MultSE was estimated. The optimal sampling range is presented independently for each locality-stratum as a light gray box. The overlapping of the boxes (dark gray) indicates the coincidence of the optimal sampling ranges between all the locations.

To characterize the hydrocarbon concentration with a precision of 0.2, the range of samples varied by locality and according to the type of hydrocarbon. The greatest effort is required for PAHs (Table 3), the range of sampling effort ranges from 4 to 14 samples per location. The localities that require more effort are Dzilam, El Cuyo, and Sisal. On the other hand, aliphatic hydrocarbons require less effort, the range is from 2 to 8 samples (Table 3). The site that requires the greatest sampling effort is El Cuyo, followed by Progreso, Sisal, and Celestún. Based on these ranges (see bottom of Table 3), we recommend sampling between 8 and 10 samples for PAHs and 6–8 samples for aliphatic hydrocarbons. As for macrofauna, these recommendations will be analyzed in terms of cost in the discussion section.

TABLE 3

Table 3. Sampling effort for estimation of hydrocarbons in sediments from intertidal strata of sandy beaches in Yucatán State, with a precision of 0.2.

Discussion

Biological Patterns

The patterns of distribution and abundance of species in sandy beaches of Yucatán State showed the following features: (1) there was a base community structure along the entire coast dominated by four species: the polychaete Polyophthalmus pictus, the oligochaete Tubificoides diazi, the isopods Excirolana mayana and E. braziliensis. Beyond these species, the community structure was highly variable in just few meters apart at the same tidal strata, but also, the effects of the tidal strata on the community structure varied within few kilometers (i.e., among sites of the same locality). (2) In general, the density of species of macrofauna was relatively low (averages between 2 and 4 species/0.01 m³) and the density of individuals was high (averages between 20 and 200 ind./0.01 m³); although both metrics were particularly higher in Dzilam. (3) The beta diversity was higher between localities with good environmental health (i.e., Dzilam and El Cuyo) but was lower between zones with regular and poor environmental health. (4) The estimated alpha diversities did not show a pattern regarding the health of the coast, they were high (>18 species) in four of six beaches. Lastly, the overall number of species reported, as well as the extrapolated, suggest that gamma diversity of the macrofauna in sandy beaches of Yucatán is within the highest known for sandy beaches.

In general, the available information about sandy beaches biodiversity around the world show a consistent pattern: richness decreases from tropical to temperate regions and from macrotidal dissipative beaches to microtidal reflective beaches (Barboza and Defeo, 2015; Fanini et al., 2020). However, besides being tropical, beaches in Yucatán are microtidal, sea breeze dominated, with a medium-grained sedimentary matrix of low content of organic matter. Therefore, considering the known patterns, the species richness reported for this region represents an atypical diversity value that motivates several ecological questions. For example, is the availability of organic matter the primary environmental filter in sandy beaches? How does the availability of organic matter interact with the morphodynamics of the beach in the stability of the macrofauna? Although these questions have been addressed previously (Defeo and McLachlan, 2005, 2013; Bozzeda et al., 2016), the pattern of diversity reported in Yucatán suggests that there is still room to address these questions in the context of beach ecology.

In Yucatán, the patterns of diversity of species do not show any geographical gradient, suggesting that local environmental/geomorphological properties might be driving the actual variation in the distribution of species. Within this context, some hypotheses that could explain these spatial patterns of species diversity are: (i) Massive seaweed accumulations as wrack during the influence of cold fronts and in dry season. The accumulation of wrack in the intertidal zone is spatially heterogeneous but widespread in this region (as evidenced in pictures of Figure 1); these accumulations in low quantities can be a source of energy and nutrients, promoting circumstantial high diversity of macrofauna (Rodil et al., 2008; Barreiro et al., 2013; Orr et al., 2014), but in large quantities could cause anoxia and large amounts of toxic leachates with negative effects on the macrofaunal community (Quillien et al., 2015). Although this process is not as evident as that of Sargassum in the Caribbean (Rodríguez-Martínez et al., 2019), it is a process that occurs annually in the coast of Yucatán State. (ii) Spatially differentiated groundwater discharges, with greater effect at the extremes of the cenote ring (Dzilam/el Cuyo and Sisal/Celestún) (Perry et al., 1995) could also generate environmental conditions for the macrofauna of the beaches. (iii) The morphodynamic conditions of each site (e.g., beach slope, swash length, erosion) might also influence in the microhabitat properties and act as an environmental filter (McLachlan and Dorvlo, 2005; Barboza et al., 2012; Barboza and Defeo, 2015). (iv) Actual anthropogenic disturbances could be reducing the beta diversity within and among some localities, especially in areas of poor and regular environmental health. To measure the plausibility of these hypotheses, it will be necessary to measure the patterns of variation in the diversity of species in times/localities with and without macrophytes wrack, with adequate spatial and temporal replication, as well as to measure environmental properties of the interstitial water (e.g., organic matter content, ammonium, redox potential, dissolved oxygen, etc.) and the morphodynamic characteristics of each site.

Regarding species distribution, at least three of the four species with wide presence along the coast of the State of Yucatán have been reported for inhabiting other sandy beaches around the world. The polychaete Polyophthalmus pictus is spread along the tropical coasts of America, but also in the Mediterranean Sea, the European Atlantic coast, South Africa, and some regions of the Pacific (Read and Fauchald, 2020). This species is recognized as a selective deposit feeder (Faulwetter et al., 2014) but also for being sensitive to organic enrichment and reported as present under unpolluted conditions (Borja et al., 2000). The oligochaete Tubificoides diazi has tanylobous prostomium and widely paired lumbricine arrangement of chaetae, this organism feeds on dissolved and particulate organic material and detritus in the surf zone. The isopod Excirolana mayana, one of the most abundant species in the present study, has been reported in the Atlantic from the Gulf of Mexico to South America, as well as in the Pacific coast of America, from Mexico to Chile (Dexter, 1976; Defeo et al., 1997; Dominici-Arosemena and Garcés, 2000). This species has been described as a very voracious predator (Brusca et al., 1995), and is one of the largest isopods in the intertidal zone (Dominici-Arosemena and Garcés, 2000). Similarly, Excirolana braziliensis has been also documented as abundant species in sandy beaches of both coasts of America as well as an important predator (Glynn et al., 1975). In Yucatán, both Excirolana species seems to compete strongly for habitat, although not explored deeply, preliminary analyzes indicate co-occurrence of only 23% on their incidences, being E. mayana the most frequent and abundant in four of six localities. Isopods of the genus Excirolana are dominant inhabitants of the supralittoral and intertidal zones of sandy beaches around the world (Omar et al., 1992). It has already been indicated that competition may be a process that promotes the exclusion and low co-occurrence of species of the genus Excilorana on sandy beaches (Defeo et al., 1997).

The study of macrofauna of beaches in Mexico dates back to the 1970s (Dexter, 1976). This first study highlights Gulf of Mexico beaches for being considerably poorer in species richness than the Caribbean beaches of Costa Rica and Panama. This appreciation changed with Méndez et al. (1985), who reported 28 species in a very extensive study (29 localities) along the >650 km of coastline of Veracruz. In the same region, Hidalgo (2017) described seasonal variability of the macrofauna, reporting 37 species. It is noticeable that the dominant species of the macrofauna in beaches of Yucatán, as well as the main assemblages’ patterns, differ from beaches of the Veracruz coast (Hidalgo, 2017). Local-scale geomorphologic and morphodynamical conditions seem to be important drivers of such differences, mainly the influence of salinity, sediment composition, wave energy, and organic input (Hidalgo et al., 2016).

Beyond these studies, and others more localized about crustaceans (Rocha-Ramírez et al., 2010; Paz-Ríos and Ardisson, 2013; Paz-Ríos et al., 2013; Rocha-Ramírez et al., 2016; Ortigosa et al., 2018), there is no information about the diversity of sandy beaches in the Gulf of Mexico. In the Caribbean, the most extensive studies are from Cuba (Ocaña et al., 2012, 2020), in which 30 species were reported; as well as some localities in Costa Rica and Panamá, with diversity varying between 3 and 15 species (Dexter, 1974). Besides these, data about the diversity of the beaches of the tropical western Atlantic is scarce. Indeed, the most important reviews about global patterns of diversity on this habitat have not included any beach of the Gulf of Mexico, neither the Caribbean (Defeo and McLachlan, 2005, 2013; Barboza and Defeo, 2015). In any case, making comparative analyzes of species richness between studies represents a challenge due to the differences in sampling effort and the extension of each study area (Nel et al., 2014). Much progress has been made using ecoregions (defined by Spalding et al., 2007) as a geographical spatial unit of comparison. For example, the richness of species in South American ecoregions ranges from 1 to 37, the peaks are located at the tropical North Brazil Shelf Province and in Guayaquil and Panama Bight (Jaramillo, 1994; Defeo et al., 2017). Other highly diverse ecoregions are from the Arabian Sea, also with richness around 21–33 species (Defeo and McLachlan, 2013; Barboza and Defeo, 2015). Therefore, from a global perspective, the species richness reported here for the coast of the State of Yucatán stands out among the highest reported worldwide.

Pollutants Patterns

Levels of hydrocarbons detected in this study are very low compared to the concentrations reported for other coastal areas of the Gulf of Mexico (GoM), both before and after Deepwater Horizon (DWH) oil spill in April 2010 (Botello et al., 2015; Bociu et al., 2019). Adhikari et al. (2016) report concentrations of ΣPAH₄₃ ranging from 68 to 160 ng/g in surface sediments collected in coastal stations in 2011 and 2013 (1 and 3 years after DWH oil spill) in the northern zone (NgoM), an area influenced by river discharges, natural oil seeps, and industrial activities including petroleum exploration/transportation (Joye, 2016). Regarding the presence of hydrocarbons in the southern part of the GoM (SGoM), De Jesús-Navarrete (1993) found concentrations of 0.8–22.6 and 34.7–79.6 μg/g of aliphatics and total PAHs, respectively, in sediments of the Campeche Sound, an area impacted by oil industry activities; Quintanilla-Mena et al. (2020) analyzed PAH content in sediment samples collected in the coasts of Mexican GoM from Tamaulipas to Yucatán, detecting the highest PAH levels in Tamaulipas (3.33 ± 5.66 μg/g) and lowest in Yucatán shores (0.04 ± 0.06 μg/g).

On the other hand, pollution studies on the coasts of the State of Yucatán are scarce. Valenzuela et al. (2005) refer maximum levels of 6.8 μg/g of aliphatics and 55.5 μg/g of total PAHs in sediments from Chelem and Progreso; Kuk-Dzul et al. (2012) report 0–3.55 μg/g of aliphatics and 0.5–5.3 μg/g of PAHs in sediments collected in four Yucatán coastal lagoons (Celestún, Chelem, Dzilam, and Ria Lagartos) in 2005; while in sediment cores collected in the sheltered port of Sisal in 2009, mean aliphatic concentrations of 1.4–9.7 μg/g and ΣPAH₁₆ of 9.9–42.8 ng/g were detected (Noreña-Barroso et al., 2017). It can be seen that the concentrations previously reported are higher than those detected in this work, including areas with the highest presence of hydrocarbons such as Progreso, Telchac and Sisal; however, it is important to consider that the sediments analyzed in those studies were collected in the coastal lagoons, which have a dynamic that makes them more vulnerable to contamination. The low levels of hydrocarbons detected on Yucatán’s sandy beaches are not surprising since the presence of these organic compounds is influenced by the potential sources of contamination, grain size (higher concentrations in fine particles) and the type of organic matter associated with the sediments (Wang et al., 2014). Sediments analyzed in this study were predominantly sandy (medium to coarse and very coarse sands) and samples were collected in areas with little industrial influence and without the presence of surface rivers. PAH concentrations in this study are several orders of magnitude lower than those considered as lower-threshold values (TELs, 312 ng/g for low MW PAHs, 655 ng/g for high MW PAHs and 1,684 ng/g for total PAHs) for marine sediment according to the quality guidelines established by NOAA’s Screening Quick Reference Tables (SquiRTs), suggesting that PAH content has a low probability of being toxic and exert a negative effect upon the benthos (Long et al., 1995; Buchman, 2008).

Beyond the oil industry’s threats, the pressure on macrofauna diversity in Yucatán could also be associated with erosion and the local strategies implemented for sediment retention. The beach erosion is a critical process in several areas of the Yucatán shoreline (Appendini et al., 2012; Cuevas Jiménez et al., 2016). In general, the transport of sediment is westward all the year; this makes the beach very vulnerable to littoral barriers, frequently used in urbanized zones to increase the width of the beach (Medellín et al., 2015; Ruiz-Martínez et al., 2016). Some of these strategies include groins, geotextile tubes, breakwaters, and beach nourishment (Ruiz-Martínez et al., 2016); however, most of these engineering strategies currently used generate down-drift erosion problems (Medellín et al., 2015; Ruiz-Martínez et al., 2016). The effects of these sedimentary dynamics on the diversity and community structure of the macrofauna have not been evaluated in this region, but it has been demonstrated in other latitudes (Walker et al., 2008; Schlacher et al., 2012; Munari, 2013). In addition, there are other sources of pollution related to the karstic characteristics of the Yucatan Peninsula, that allow pollutants generated by activities carried out inland to permeate, flow and enter to the coast by submerged groundwater discharges; some evidence of that is the presence of human fecal material (coprostanol and epicoprostanol) in sediments surrounding the springs located in Dzilam de Bravo (Kantun-Manzano et al., 2018), as well as the results reported by Arcega-Cabrera et al. (2015) concerning heavy metal pollution in Chelem Lagoon connected to inland anthropogenic activities along with local factors.

Sampling Effort

For macrofauna, we recommend that the sampling effort for future ecological assessments on sandy beaches of Yucatán should be of eight cores per stratum at each site, with four sites (separated by 2–3 km) for each one of the six localities, at least three times annually, considering the temporal variation of beach morphodynamics (i.e., erosion, accretion and accumulation of algae) (Medellín and Torres-Freyermuth, 2019). This full design would improve the precision with which the species composition is estimated along the coast of Yucatán, as well as the natural temporal variability. It also will allow to obtain a total sampling area for each site of about 3.6 m²; a bit above the minimum recommended by Schlacher et al. (2008) for sites in microtidal beaches with less than 50 m beach width (i.e., between 2 and 3 m²). Especially, this full design would allow bBACI analyses (beyond Before – After × Control – Impact sensu Underwood, 1992) in the context of oil spills or any other major disturbance that might occur on the beaches of Yucatán (Bejarano and Michel, 2016). However, supporting this sampling effort for the six beaches might be not affordable in terms of the time and personnel required to process so many samples. In such a situation, three strategies could be applied to reduce the sampling effort: (1) reclassifying littoral zones into dry zone and swash zone, as recommended by MBON P2P (Pole, 2019) for microtidal beaches; (2) dispensing to sample in two or three localities with proved redundant information (e.g., Celestún, Progreso, Telchac); and (3) reducing to two annual samplings. Hence, keeping the most biodiverse localities (i.e., Dzilam and El Cuyo) and one of the remaining locations in the sampling design (e.g., Sisal) would ensure that the spatial patterns of biodiversity will be represented and measured twice a year, although potential reference localities would be sacrificed, as well as the temporal variability, which would imply a weakening of the bBACI analysis. Put in numbers, the full design would involve taking and processing 576 samples each time (1,728 samples per year), while the reduced sampling would involve 192 samples (384 samples per year). As a reference, this study involved processing and analysis of 225 sediment samples.

For pollutants, statistical estimations indicated sample sizes of 8–10 samples for PAHs and 6–8 samples for aliphatic hydrocarbons at each stratum. Based on these, we recommend taking eight sand samples at each stratum for both kinds of pollutants, for each of the four sites of each of six localities. As for macrofauna, this sampling effort generates several benefits. The first is the high statistical precision (and power) to detect changes in a highly variable system, using a bBACI analysis. The second is the feasibility of building regression models between the biological properties of the macrofauna concerning the variations of the pollutants in sediments. Although this is not informative in the actual circumstances (due to the extremely low concentrations of pollutants), it will be highly informative in case of an oil spill in the future. The great limitation of this sampling design is the economic feasibility; estimating PAH and aliphatic hydrocarbons is expensive, in orders of magnitude higher than the processing of biological samples. Hence, a reduction of the sampling effort as the one proposed for the macrofauna could also be applied for contaminants, with the disadvantage of losing information on potential reference locations.

Conclusions

In this work, we demonstrate that the diversity of macrofauna is heterogeneous along the coast but especially high on the east side of the state. We also demonstrate that the current levels of PAHs and aliphatic hydrocarbons in beaches are low. Both situations indicate that it is an appropriate time to start long-term monitoring of the environmental health of Yucatán beaches, even more now when energy development (mainly oil & gas) in the Yucatán platform is imminent. It is necessary to advance in the understanding of the temporal patterns of diversity and abundance of species to complement the ecological baseline necessary to identify environmental impacts. The sampling designs we suggest are based on statistical precision and in recommendations previously formulated by experts in beach ecology. We also propose alternatives to reduce sampling effort concerning the ideal design considering possible resource limitations (e.g., funds, personnel, etc.), which, as in most cases, is usually the limiting factor for robust environmental impact assessments. The sampling protocol developed by the MBON Pole to Pole will allow comparisons to be made with other beaches throughout the Americas.

Further research directions must include the role of morphodynamics (natural and human-induced) and temporal variability in the supply of energy and nutrients as main environmental drivers of macrofauna diversity. Besides, keeping the focus of biological and pollutant spatiotemporal patterns, and developing long-term data sets with international standards, are within the expectations for modern integrative beach ecology and assessment of environmental impacts.

Data Availability Statement

The datasets presented in this study can be found simultaneously in OBIS or Zenodo data repositories with the respective accession codes: OBIS – https://obis.org/dataset/5700fdc3-956f-4e33-903b-e66ba38d980c and ZENODO – doi: 10.5281/zenodo.3771828.

Author Contributions

NS, MM, EN-B, GH, and EG-C contributed to the conception and design of the study. EG-C, GH, and RC-C did fieldwork. GH, RC-C, and MM-R processed biological samples. JQ-D and EN-B did granulometric and chemical analyses. GH, JQ-D, and EG-C organized the database. EG-C performed statistical analyses. EG-C wrote the first draft of the manuscript, with the collaboration of NS, MM-R, and EN-B. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This project was funded by the 2018 CONACyT 293354 grant “Laboratorio Nacional de Resiliencia Costera.” The continuity of this research is financed by the grant IA206320 (PAPIIT, Universidad Nacional Autónoma de México) under the supervision of EG-C. Funding was also provided by the Harte Research Institute for Gulf of Mexico Studies and the Harte Charitable Foundation through the “Biodiversity of the southern Gulf of Mexico” project under the supervision of NS and the support of the BDMY research group. NS holds the Furgason Fellowship International Chair for Coastal and Marine Studies in México.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Special thanks to the coordinators of the Pole to Pole of the Americas network, as well as to all its members for their positive comments and constructive criticism of this work in the workshop held in Puerto Morelos, Mexico, in 2018. We extend our gratitude to Dr. Frank Ocaña, who actively participated at an early stage of the design of this project. This work benefited greatly from important insights and comments given by reviewers.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.589656/full#supplementary-material

Footnotes

References

Adhikari, P. L., Maiti, K., Overton, E. B., Rosenheim, B. E., and Marx, B. D. (2016). Distributions and accumulation rates of polycyclic aromatic hydrocarbons in the northern Gulf of Mexico sediments. Environ. Pollut. 212, 413–423. doi: 10.1016/j.envpol.2016.01.064

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, M. J., Crist, T. O., Chase, J. M., Vellend, M., Inouye, B. D., Freestone, A. L., et al. (2010). Navigating the multiple meanings of (diversity: a roadmap for the practicing ecologist. Ecol. Lett. 14, 19–28. doi: 10.1111/j.1461-0248.2010.01552.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, M. J. (ed.). (2017). “Permutational multivariate analysis of variance (PERMANOVA),” in Wiley StatsRef: Statistics Reference Online, (Hoboken, NJ: John Wiley & Sons, Ltd). doi: 10.1002/9781118445112.stat07841