Testing marine regional-scale hypotheses along the Yucatan continental shelf using soft-bottom macrofauna

Different hypotheses related to the regional-scale configuration of the Yucatan Continental Shelf (YCS) between the Gulf of Mexico (GoM) and the Caribbean Sea have been proposed. Hypotheses regarding its regional boundaries include: (i) an ecoregional boundary at Catoche Cape, dividing the Western Caribbean and the Southern GoM ecoregions; and (ii) a boundary within the Southern GoM ecoregion at 89°W, separating the West and Mid-Yucatan areas. We tested the hypothesis of no variation in benthic macrofaunal assemblages between regions delimited by the former boundaries using the species and functional traits of soft-bottom macrofauna. We considered that the depth and temporal environmental dynamics might interact with regional variations, generating complex benthic community patterns. The data were collected over five years (2010–2012, 2015–2016) at 86 stations (N = 1, 017 samples, 10–270 m depth), comprising 1,327 species with 45 combinations of functional traits. The variation in species composition and functional trait assemblages were both consistent with the occurrence of three separate regions in the Yucatan Peninsula (West Yucatan, Mid-Yucatan and Western Caribbean). This regional configuration was consistent with changes in assemblage structure and depth zonation as well as temporal variation. Along with spatial and temporal variation, diversity diminished with depth and different regions exhibited contrasting patterns in this regard. Our results suggest that the spatial and temporal variation of soft-bottom macrofauna at YCS demonstrate the complex organization of a carbonate shelf encompassing different regions, which may represent transitional regions between the Caribbean and the GoM.


INTRODUCTION
Ecoregions are defined as areas with relatively homogeneous communities that are distinct from adjacent systems (Spalding et al., 2007) and that are also affected by environmental conditions, such as currents, upwellings, primary productivity and sediment, among others. The delimitation of ecoregions enables us to understand biogeographical processes, large-scale diversity patterns and their environmental drivers (Paganelli, Marchini & Occhipinti-Ambrogi, 2012;Williams et al., 2015) and to define conservation priorities (Selig et al., 2014;Tear et al., 2014;Beger et al., 2015).
A comprehensive definition of the configuration of marine ecoregions has been proposed by Spalding et al. (2007) based on biogeographical data and environmental conditions. Since then, new evidence has emerged to either support or challenge this means of accounting for changes in the composition of assemblages (Williams et al., 2015) and even benthic deep-water assemblages (Hernández-Ávila, 2014;Hernández-Ávila et al., 2018) as well as genetic connectivity (DeBoer et al., 2014). Moreover, one current debate focuses on the delimitation of some ecoregional boundaries (Carpenter et al., 2011). For instance, a distinct regional configuration has been proposed for the Gulf of Mexico based on chlorophyll-a concentration patterns (Salmerón-García et al., 2011), which may help to redefine ecoregions in smaller areas. Nevertheless, the effects of depth on spatial variation at a regional scale remain poorly explored. There is a broad consensus that the environmental drivers that influence the distribution of benthos change with depth, leading to a gradient effect on the dynamics of the abundance and distribution of benthic assemblages from shallow to deep water (Zajac, 2008;McArthur et al., 2010). Testing different regional variation hypotheses and including depth in spatial analyses of marine benthos may assist in the identification of ecological variation at regional scales.
Soft-bottom marine assemblages are useful models for testing meso-and large-scale variation hypotheses because they can respond to changes in the main environmental factors used to define marine biogeographical areas, such as oceanographic conditions, depth, organic carbon inputs and sedimentary provinces (Snelgrove & Butman, 1994). Ecosystem functioning is associated with the biodiversity of soft-bottom habitats (Snelgrove et al., 2014;Strong et al., 2015), which can be studied according to the diversity of species and variations in their functional biological traits (Menezes, Baird & Soares, 2010). Analyzing biological traits may prove helpful in detecting ecosystem patterns and functioning and anthropogenic disturbances (Bremner, 2008;Suding et al., 2008;Vinagre et al., 2017) as well as the response of an assemblage to environmental conditions (Menezes, Baird & Soares, 2010;Paganelli, Marchini & Occhipinti-Ambrogi, 2012; Van der Linden et al., 2012).
According to Spalding et al. (2007), the Yucatan Shelf separates two marine ecoregions: the Southern Gulf of Mexico (GoM) and the Western Caribbean. The boundary between these two regions follows the northeastern contour of the Yucatan Shelf from Catoche Cape (21.6 • N 87.1 • W), which is bounded in the north by the Greater Antilles ecoregion. The Southern GoM ecoregion extends up to the Tamaulipas coast, where it meets the Northern GoM ecoregion. However, Salmerón-García et al. (2011) have proposed a method of regionalizing the GoM based on the regional pattern of chlorophyll-a, which suggests that both the Northern and Southern GoM ecoregions can be subdivided into smaller regions. The regional setting of Salmerón- García et al. (2011) includes the separation of coastal Yucatan areas (<approximately 50 m deep) between the Tabasco-Campeche shelves and a section of the Campeche Bank extending up to approximately 89 • W. Although different regional settings of the Yucatan area have been accepted in biogeographical analyses of marine assemblages (Francisco-Ramos & Arias-González, 2013;Williams et al., 2015), a degree of uncertainty exists regarding the marine communities along the Yucatan Peninsula. The aim of this paper is to test the hypothesis of no variation in macrofaunal assemblages between regions at the Yucatan Continental Shelf (YCS) through analyzing their patterns in terms of species composition and functional traits. The combined analysis of species and functional assemblages may offer a more integrated view of large-scale variations in marine communities by identifying different variation patterns between ecoregions with both depth and time. Given that the extent of coastal influence decreases with depth (Weissberger et al., 2008;Zalmon et al., 2013), assemblage structures may differ among ecoregions due to bottom topography. The samples for statistical analysis comprise multiannual data (2010, 2011, 2012, 2015 and 2016) collected along the depth range of the YCS that consider both temporal variations and the effects of depth on the composition of assemblages.

Study area and regional hypotheses
The YCS is a carbonate shelf extending over 100-300 km from the shoreline, covering an area of approximately 57,000 km 2 and composed of sedimentary deposits of mediumto fine-grained carbonate sand (Balsam & Beeson, 2003;Logan et al., 1969). Owing to the karstic nature of the Yucatan Peninsula, most of the freshwater input (8. 6×10 6 m 3 km −1 yr −1 ) in the area comes from groundwater (Hanshaw & Back, 1980), affecting nutrient input and phytoplankton communities (Slomp & Van Cappellen, 2004;Troccoli-Ghinaglia et al., 2010). The Yucatan Shelf hosts various coral reefs near the coast and at the shelf margin (Zarco-Perello et al., 2013) as well as coastal lagoons and mangroves (Hernández-Guevara, Pech & Ardisson, 2008) along the shoreline. The major currents in the area include the GoM loop current, which primarily affects the eastern and northern sections of the Yucatan Peninsula, and westerly currents along the west coast (Zavala-Hidalgo, Morey & O'Brien, 2003). Historically, the YCS has been subjected to hurricanes (Boose et al., 2003).
The regional-scale differences tested in this study include the division proposed by Spalding et al. (2007) between the Western Caribbean and the Southern GoM (Fig. 1). The northwestern boundary of the Southern GoM ecoregion is located at 23 • N 97.8 • W and is characterized by lower variation in the sea surface temperature during the winter than that observed in the Northern GoM ecoregion (Spalding et al., 2007;Wilkinson et al., 2009). The Southern GoM ecoregion also includes a subdivision between the YCS and Tabasco Shelf at around 92 • W, in the Campeche Canyon (Lara-Lara et al., 2008;Wilkinson et al., 2009). According to this configuration, the western section of the YCS from Catoche Cape is considered as a single regional unit within the Southern GoM ecoregion (Yáñez Arancibia & Day, 2004;Zavala-Hidalgo & Fernandez-Eguiarte, 2007;Lara-Lara et al., 2008;Wilkinson et al., 2009). An additional regional division of the Southern GoM, based on spatial variations in chlorophyll-a concentration, has been proposed by Salmerón-García et al. (2011) and includes a boundary at around 89 • W, thereby separating a region in the centre of the YCS from the western section that includes Campeche Bank and the Tabasco Shelf. However, this regional setting does not include the boundary at Catoche Cape proposed by Spalding et al. (2007). Moreover, Salmerón-García et al. (2011) have proposed another ecoregion in the outer section of the YCS, defined by processes associated with the transition between coastal and oceanic zones (Ruíz-Castillo et al., 2016). We compared macrofaunal assemblages across three regions proposed in previous works: (a) West Yucatan, as part of the Southern GoM region proposed by Spalding et al. (2007), including locations west of 89 • W according to Salmerón-García et al. (2011); (b) Mid-Yucatan, located from the eastern boundary of the West Yucatan region to the Catoche Cape boundary, suggested by Spalding et al. (2007); and (c) the Western Caribbean, located east of Catoche Cape. By comparing regional differences in soft-bottom macrofaunal assemblages, we tested the hypothesis of a lack of regional-scale variation in benthic macrofaunal assemblages on the Yucatan Peninsula according to the boundaries proposed by Spalding et al. (2007), Salmerón-García et al. (2011) or both. All samples were collected using a 0.25-m 2 Smith Mcintyre sediment grab, collecting three 10 cm random cores per grab at each station. Benthic macrofauna were recovered from the first 15 cm of the sediments using a 500 µm sieve. Once recovered, the organisms were pre-fixed using 15% magnesium chloride solution for 15 min and then fixed in a 4% formaldehyde buffer solution. Benthic macrofauna were identified to the lowest taxonomic level possible, typically to species or morphospecies level. All procedures from data sampling to database analysis, including taxonomic determination, were conducted following a quality assurance/quality control (QA/QC) protocol (D. Pech, 2018, unpublished data) developed by the Biodiversidad Marina y Cambio Climático (BIOMARCCA) laboratory.

Identification of functional traits
Five functional traits commonly used for benthic fauna (Costello et al., 2015) were selected to represent the ecological functioning or the response of taxa to their environment. These traits were divided into categories (Table 1) and used to generate a binary taxa × trait matrix, where one (1) represented the fit of the taxa to the dominant trait category and zero (0) represented ''no'' fit (Beauchard et al., 2017). Traits were mostly identified based on a literature review conducted for each species or close relative taxa (genera or within family levels). Some traits were estimated based on direct observations (i.e., the presence and type of skeleton or feeding structures) and by consulting experts. One category per trait was selected to identify species trait combinations, with the exception of feeding trait, which presented one or two categories (e.g., predator and scavenger). Species with records of more than two feeding categories were included as omnivorous. The combination of trait categories was used to construct a trait combination × sample matrix by merging the abundance of species fitted in the same combination of functional traits for each sample.

Data treatment Estimation of taxa and functional trait contributions to the assemblage
To combine the contributions of the taxa to the general community composition based on their presence, abundance and richness, we calculated a taxon contribution index (TCI) as follows: where pr i , ar i , and rr i were the corresponding presence, abundance and richness ratios of ''i'' taxa, respectively, estimated as: rr i = species i species (4) The TCI i values ranged from 0 to 1. The maximum score was obtained when taxon ''i'' was always present in the sampling data set and contained all possible species and individuals. Taxa with higher species and levels of abundance could exhibit high TCI scores, but this was also true of taxa with few species where their occurrence in the sampling set was high. The TCI was not estimated in the place of traditional indices, but rather as a general combined reference for the taxon's relevance within the community composition. The index was estimated using Class-level taxonomic references, with the exception of crustaceans, which we included at the Order level. Similarly, to estimate the relative relevance of trait combinations, the functional trait index (FTI) was estimated based on the ratio of species richness and the presence and abundance of each combination of functional traits. The estimations were based on Eqs.

Estimation of species diversity
To compare species richness between depth ranges, accumulation species curves (obtained by sampling and the Chao 1 estimator) were plotted as functions of depth ranges using ESTIMATES (Colwell, 2013). Estimations were conducted per depth interval (<50 m; 50-100 m; >100 m), extrapolated to 400 samples (for the species cumulative curve) and based on 100 permutations. Extrapolations and confidence intervals were estimated using the derivation of the Bernoulli product model (Chao et al., 2009;Clowell et al., 2012). Estimations were conducted using data obtained from all cruises. Differences in species richness as a function of depth were tested by conducting t-tests based on the parameters obtained by ESTIMATES. In addition, we compared the confidence intervals of each curve (Colwell, Mao & Chang, 2004;Clowell et al., 2012).

Testing regional-scale variations in species composition and functional traits
Two similarity matrices were calculated using the Bray-Curtis Index based on both the abundance of each species and the abundance of each combination of functional traits. The null hypothesis of non-differences in the assemblage compositions between regions, depth zones and years was tested by conducting permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001) for both cases. The PERMANOVA tests included Region and Year as factors, with sampling depth (10-270 m) as the covariable. Region (Western Caribbean, Mid-Yucatan, West Yucatan) and Year (2010Year ( , 2011Year ( , 2012Year ( , 2015Year ( , 2016 were included as fixed factors. In addition, sampling station was included as a nested factor (16-41 stations per Region). Three pseudoreplicates were considered for each sampling point. The null model was tested using 9,999 permutations of residuals under the reduced model. For fixed factors, PERMANOVA t-tests were used to conduct pairwise comparisons. The effect of depth on assemblage structure was also explored using canonical analysis of principal coordinates (CAP), using station centroids against mean depth. The significance of the first-square canonical correlations was tested using 999 permutations (Clarke & Gorley, 2015). The interactions between temporal variation (Year) and regional variation (Region) in terms of assemblage composition were explored using second-stage analyses (Clarke et al., 2006) for both species and functional trait assemblages. Regarding the absence of an interaction effect between Year and Region, high second-stage correlations between regional matrices of temporal variation were expected (Clarke et al., 2006;Clarke & Gorley, 2015). The distribution trends of species assemblages and functional trait assemblages were represented by the bootstrap average of centroids in multidimensional scaling plots (MDS) using 95% confidence intervals and 30 bootstraps per group (Clarke & Gorley, 2015).

Identifying influential species
The original species assemblage was reduced to a subset of species that could represent the overall pattern using a BVStep routine (Clarke & Warwick, 1998;Clarke & Gorley, 2015). Similarity data matrices with different combinations of species subsets were compared with the similarity matrix containing all taxa to identify the smallest group of species that could best describe most of the pattern in the full data set. An iterative species exclusion-inclusion process was performed during the BVStep routine to obtain the minimum number of species from the species subset that could represent a similarity matrix that was highly correlated with the original similarity matrix. The Spearman's rank coefficient was used to compare similarity matrices using a threshold of ρ > 0.90 as the minimum level of correlation between all taxa and species subset matrices. BVStep analyses were performed several times (over 30 trials) to ensure that the best species subset was selected.
The combination of trait categories produced 45 functional trait combinations for the YCS soft-bottom macrofauna. Five functional trait combinations accounted for most of the specimens in terms of species and occurrence ( Fig. 2B): (i) organisms with a soft, flexible body and limited motility, a deposit feeding mode, providing no parental care for their eggs and producing planktonic larvae; (ii) organisms similar to the first group, but with a predator feeding mode; (iii) organisms with a fragile exoskeleton and relatively high motility, a deposit feeding mode and direct development, with limited dispersion of their propagules; (iv) organisms with soft, flexible bodies and limited motility, with planktonic larvae and omnivorous feeding; and (v) shelled organisms with limited motility, a filter feeding mode and planktonic larvae.

Trends in species richness
Overall, the cumulative species diversity (1,526.2 ± 24.5 spp., extrapolated to 1,500 samples) was slightly lower than the Chao 1 estimator (1,742.7 ± 52.5 spp.), with both exhibiting an almost asymptotic pattern. Higher species richness was estimated at shallower depths Colors in trait codification refers to structural fragility (blue), motility (green), living position (brow), reproduction/development mode (red) and feeding (black). Reference codes for trait categories in Table 1.
The decrease in diversity with depth at the YCS was also consistent across all sampling years.

Changes in species composition and functional traits along the Yucatan continental shelf
According to the PERMANOVA test, variations in assemblage structure based on species composition and functional traits exhibited consistent responses to spatial (region, depth) and temporal (year) factors and their combinations ( Table 2). The significant correlation between station centroids and depth (species assemblage: 2 = 0.832; functional trait assemblage: 2 = 0.664; p = 0.001 both cases) suggest an effect of depth on species and functional assemblages. In both types of assemblage, the effect of depth differed by Region and Year, denoting varied zonation patterns. Significant interactions were detected between region and year. A posteriori analyses indicated significant changes in species composition between all years in the three ecoregions. Functional trait assemblages also presented variations in most of the pairwise analyses between years, with different patterns between regions (Table 3). The two-stage analyses of region per depth range centroids exhibited low correlations between years (Spearman's correlation, species assemblage ρ < 0.428, functional trait assemblage ρ < 0.518), supporting the occurrence of different temporal variation patterns occurring at each region and depth range. The variation of the temporal trends between the assemblages of the regions constituted the interactive effects detected between Region and Year. The a posteriori analysis performed for each depth range revealed that the species assemblage showed differences between regions at all depths (Table 4). For the functional trait assemblage, significant differences were found for the upper shelf, but the spatial variation decreased at lower sections. At shallow depths (<50 m), significant differences between the Western Caribbean and both Mid-and West Yucatan were found. At middepth (50-100 m), differences between Mid-and West Yucatan were detected, while no differences were found between these regions and the Western Caribbean. Beyond 100 m, no differences were detected between the three regions ( Table 4). The MDS plot based on centroids indicated that the regional differences observed in the species and functional trait assemblages at the shallow and mid-depths diminished at depths exceeding 100 m (Fig. 4).

Variation in functional traits by ecoregion and depth range
Consistent with the PERMANOVA test, variations in the distribution of functional traits were observed. In general, reductions in the abundance of most trait combinations were observed from shallower to greater depths. The depth pattern differed between trait combinations and regions, as observed in Table 5 and Table S2. Significant variations in trait combinations between the Western Caribbean and both Mid-Yucatan and West Yucatan were detected at 10-50 m (Table 4), where the Western Caribbean was dominated by infaunal soft-bodied predators while West Yucatan was dominated by two types of infaunal deposit feeders (organisms with soft-bodies and limited motility with indirect development; and exoskeleton-bearing, motile organisms with direct development). In Mid-Yucatan, infauna were dominated by deposit feeders, predators and herbivores (including soft-bodied species with limited motility and indirect development, and motile, exoskeleton-bearing species that undergo direct development). At a depth of 50-100 m, the differences between Mid-and West Yucatan were associated with an increase in soft-bodied predators, herbivores and exoskeleton-bearing detritivores in the former.

DISCUSSION
The overall species diversity at the YCS was much higher than that of other soft-bottom macrofaunal shelf communities (Table 6). Only three locations with a comparable scale of diversity (>1,000 species) were identified: the Campos Basin of Rio de Janeiro (Zalmon et al., 2013), the Monterey Bay Continental Shelf (Oliver et al., 2011) (both of which have upwelling systems, like the YCS) and the carbonate shelf at the northern limit of the GoM (SW Florida) (Phillips, David & Keith, 1990). No other locations with high carbon input from riverine or upwelling origins exhibited a similar level of diversity. In the case of the YCS, the occurrence of different types of habitats (such as coral reefs, seagrass beds, coastal lagoons and groundwater inputs) could contribute to species diversity through some Table 5 Abundance (ind 0.1m -2 , mean ± se) of 10 relevant taxa from the representative sub-set of 60 the taxa that reproduce the observed patterns, and the 10 most abundant combination of functional traits of the macrofauna at YCS. Abundance is estimated according to depth range and regions. Abbreviations in brackets, Pol, Polychaeta; Amp, Amphipoda; Tan, Tanaidacea; Cop, Copepoda; Gas, Gastropoda; Sip, Sipuncula. Colors in functional trait codification refers to structural fragility (blue), motility (green), living position (brown), reproduction/development mode (red) and feeding mode (black). Reference codes for functional traits are the same used in Table 1. Information on the subset of the 60 representative taxa are show in Table S1 and from the 20 most abundant combinations of functional traits in Table S2.  26.2 ± 6.7 5.9 ± 1.4 4.9 ± 1.5 Plesiolembos ovalipes (Amp: Aoridae) 9.26 ± 5.7 1.9 ± 1.2 2.8 ± 2.0  degree of connectivity between habitats. Most species in the three different assemblages observed here (69%) undergo indirect development with larval dispersal. The Yucatan Current could enhance larval connectivity by facilitating the recruitment of many species. The general structure of the soft-bottom macrofauna was dominated by polychaetes and peracarid crustaceans, followed by molluscans (Bivalvia and Gastropoda), decapod crustaceans and ostracods, similar to the general patterns of the benthic communities in continental shelf habitats, although the lack of species-level classification conducted in previous assessments near our study area (Hernández-Arana et al., 2003;Escobar-Briones & Jimenez-Guadarrama, 2010) precludes direct comparison. Similar patterns have been detected using family-level information in previous studies on the western section of the YCS, with common dominant taxa including polychaetes from the families Spionidae, Syllidae, Paraonidae, Aspidosiphonidae, Sabellidae and Lumbrineridae (Hernández-Arana et al., 2003).
Consistent functional trait combinations that significantly contributed to the assemblages include: (i) infaunal organisms with a soft-flexible body, limited motility, larval dispersion and diverse feeding modes, such as deposit feeders, scavengers and predators (mainly polychaetes); and (ii) infaunal, exoskeleton-bearing motile species that undergo direct development with the deposit feeding or scavenging trophic mode (mainly peracarid crustaceans). Furthermore, the number of functional trait combinations decreased with depth, similar to the species diversity patterns.
The results for both the species and the functional trait assemblages were consistent with the regionalization of soft-bottom macrofauna at YCS. Biogeographical boundaries are related with shared species and diversity, but recent studies have highlighted the relevance of using a trait-based approach to describe marine biogeography (Brun, Payne & Kiørboe, 2016). Here, we demonstrate that trait analysis is a useful surrogate for depicting regional-scale differences in benthic assemblages, particularly for shallow depths. Despite detecting regional variations in species assemblages at all depth ranges, macrofaunal assemblages based on functional traits only exhibited regional variations in the shallow section of the shelf (10-50 m). This may have been a result of a regional variation in species with the same combination of functional traits. By including more traits and categories, it may be possible to achieve a greater resolution of functional trait combinations, facilitating the detection of spatial and temporal variation with higher accuracy.
The regional configuration of the soft-bottom benthic community represented a combination of the previous regional boundaries proposed for the area. The boundary proposed at Catoche Cape by Spalding et al. (2007), Lara-Lara et al. (2008) and Wilkinson et al. (2009), separating the Western Caribbean ecoregion, was consistent with our analyses. However, our data also suggest that the western section of YCS from Catoche Cape is not a single regional unit. Variations between Mid-and West Yucatan were identified, supporting these regions' separation according to Salmerón-García et al. (2011). This result suggests that considering the Southern GoM ecoregion as a single unit may underestimate the ecological variability of the GoM. Detecting regional-scale variations in marine assemblages associated with environmental drivers in the Southern GoM would contribute to redefining smaller ecoregions within the natural regional-scale variation of the GoM. In the current case of the YCS, each area has a particular range of temporal and spatial (depth) variability, supporting the separation of three distinct ecoregions.
The effect of depth on assemblage composition differed among regions and years. The different benthic assemblage compositions could be associated with the varied shelf configurations observed along the YCS. Indeed, there is a short shelf close to the continent and a stepped slope in the Western Caribbean section. However, a very wide shelf is also present that extends 200-300 km from the continent with an almost steady slope at the West and Mid-Yucatan regions. Environmental conditions interacting with the physical configuration of the shelf may drive the different temporal dynamics of benthic communities, explaining the changes in temporal patterns and variations with depth between ecoregions.
Comparing the composition of assemblages at each depth range, both species and functional trait assemblages exhibited variations in assemblage structure between shallow (<50 m depth), mid-(50-100 m) and >100 m depths. According to Spalding et al. (2007), the ecoregions proposed for marine areas apply to coastal and shelf environments up to a depth of 200 m. Nonetheless, current data suggest that deeper shelf regions (>50 m deep) may exhibit less regional-scale variation in functional trait assemblages than shallow waters, although variations in assemblage compositions in the area were detected between the three ecoregions when considering the total depth range of the shelf. Comprehensive studies on marine biogeography have focused on large, global-scale depth ranges (Watling et al., 2013), but conducting a more precise separation of depth intervals for shelf environments might improve ecoregional classification, as demonstrated here.
The variations in assemblage composition may have been associated with differences in the water masses of the YCS. Water masses have already been suggested as the main drivers of the biogeographic structure of benthic shallow-water and continental-shelf species worldwide (Bellanger et al., 2012). The thermohaline conditions of the central coast of the Yucatan Peninsula may define major spatial oceanographic variations in areas where two types of bottom-water masses occur at shallow depths, including i) the local water mass (Yucatan Sea Water or YSW, temperature of 26-31 • C and salinity of 36.4-36.8), and ii) the upwelled Caribbean Subtropical Underwater (CSUW, <23 • C, salinity 36.25-36.75) (Enriquez et al., 2013). CSUW upwelling occurs at Cape Catoche and its effects extends offshore to the west up to 89.5 • W, while the YSW occurs at both the central and western sections of the YCS, with a strong presence west of 89 • W. Bottom water with low salinity was observed in the eastern section of the YCS due to the influence of submarine groundwater discharges (SGD) from the Holbox fractures occurring between 88 and 87.5 • W (Pope, Ocampo & Duller, 1993). The combined effects of the CSUW upwelling and SGD may generate high primary productivity in the central section of the YCS, affecting the assemblages of soft-bottom communities, particularly those in the inner section of the shelf (Ruíz-Castillo et al., 2016). The westward circulation west of the 89 • W upwelling may drive the extension of the upwelling towards the inner shelf, but with less influence on benthic communities. We suggest a conceptual model based on our results (Fig. 5) that indicates that soft-bottom communities in West Yucatan may be affected by YSW at shallow depths (0-50 m), Gulf Common Water (GCW) occurring from 0 to 250 m and by the wind-driven westward circulation of upwelled water (Ruíz-Castillo et al., 2016) and local SGD.
The non-differences between assemblages in the deeper section of the shelf (>100 m) may owe to the more homogeneous offshore oceanographic conditions. The lack of differences between deep water masses may generate more homogeneous bottom conditions than those observed at shallow depths. Moreover, local upwellings typically exert important inshore effects on bottom communities, but their extension to offshore waters depends on inshore-offshore transport (Ruíz-Castillo et al., 2016). The influence of inshore upwellings on deep offshore shelf communities can be lower than that on shallow-water communities. Beyond a depth of 100 m, SGD has little effect (if present) and bottom habitats in the area are affected by CSUW (150-250 m) and GCW (0-250 m) (Enriquez et al., 2013). According to Vázquez (2000) and Enriquez et al. (2013), GCW is the CSUW affected by mixing processes and local changes.
The delimitation of the study area did not allow us to determine the western limit of the West Yucatan region. However, owing to the scale of the spatial variations observed, its western limit could be located at the limit of the YCS, close to the Tabasco Shelf. In the western section of the Yucatan Continental Shelf (near Campeche Canyon, around 92 • W), major differences in sedimentary settings have been observe, as they  (Spalding et al., 2005). Also included is a suggested western limit base on the transition of carbonate sediments on the shelf and their effects on macroinfauna composition (Hernández-Arana et al., 2005). Gray arrows represent the loop current (dashed) and dominant currents on the shelf according to Kjerfve (1994), Zavala-Hidalgo, Morey & O'Brien (2003) and Lie-Yauw (2005). On inserted squares are species from the 60spp subset that were more abundant, comparing other regions. Brackets indicate depth range: s, 10-50 m; m, 50-100 m; d, >100 m. 1 (Hernández-Arana et al., 2005), 2 (Salmerón- García et al., 2011), 3 (Spalding et al., 2007), 4 (Pope, Ocampo & Duller, 1993), 5 (Ruíz-Castillo et al., 2016), 6 (Kjerfve, 1994. Full-size DOI: 10.7717/peerj.8227/ fig-5 shift from a carbonate shelf to terrigenous sediments (Carranza-Edwards, Rosales-Hoz & Monreal-Gómez, 1993;Hernández-Arana et al., 2003;Hernández-Arana et al., 2005). These environmental conditions are also associated with shifts in the composition of soft-bottom communities (Hernández-Arana et al., 2003). This suggests that changes in the sedimentary provinces close to Campeche Canyon could generate additional biogeographical subdivisions in the western section of the YCS. There is evidence to suggest that West Yucatan macrofauna change according to the dominant climatic periods occurring in the area (Hernández-Arana et al., 2003;Hernández-Guevara, Pech & Ardisson, 2008). Here, samples were collected during the late rainy seasons of 2010, 2011 and 2016 and the beginning of the cold front periods of 2012 and 2015. Some differences detected between the cruises were probably caused by differences between years and periods (between 2011-2012 and 2015-2016). Although seasonal variations in the coastal and shallow regions of the YCS are usually associated with continental run-off, wind stress and sea surface temperature (Hernández-Arana et al., 2003;Hernández-Guevara, Pech & Ardisson, 2008;Kuk-Dzul, Gold-Bouchot & Ardisson, 2012), climatic variations generate oceanographic conditions that can affect deeper waters. The oceanographic conditions of the Southern GoM during the rainy season include low wind stress (<5 m s −1 ), low mixed layer depth (<40 m) and low net primary production (<250 mg C m −2 day −1 ), while cold front conditions include high wind stress (>6 m s −1 ), a greater mixed layer depth (>80 m) and high primary productivity (>350 mgC m −2 day −1 ) (Müller-Karger et al., 2015). These contrasting scenarios may be associated with changes in carbon inputs to the bottom shelf, affecting the benthic infauna. Disentangling the potential seasonal variation from annual variation would allow for a better understanding of the dynamics of the benthic communities on the YCS.
In addition, annual variations were detected between the samples collected during the rainy season (2010 and 2011) as well as during the early cold front period (2012 and 2015). Although there is evidence that GoM shelf macrofauna change according to seasonal variations (Posey et al., 1996;Posey et al., 1998;Hernández-Arana et al., 2003), annual variations may also occur within the same period, but this aspect has been relatively less explored. In addition to seasonal variations, interannual climatological variations in the GoM are responsible for variations in phytoplankton biomass and primary production (Müller-Karger et al., 1991;Martínez-López & Zavala-Hidalgo, 2009;Müller-Karger et al., 2015). Variations in biomass generated at the sea surface would affect the input of carbon to the benthic shelf, influencing benthic communities. The potential seasonal and interannual variations associated with oceanographic conditions and phytoplankton biomass must be addressed by coupling environmental data and oceanographic models with the structure of species and functional trait assemblages.

CONCLUSIONS
Our findings have demonstrated spatial, temporal and depth variations in the distribution of soft-bottom macrofaunal assemblages along the YCS. The spatial and temporal variations indicated the complex organization of carbonate shelf communities, which were previously believed to be relatively homogeneous environments, where major spatial and temporal changes can occur. The ecoregional boundaries of the Southern GoM (according to general environmental conditions) should be re-evaluated, considering the more precise delimitation of environmental and community assemblage variations. Our results suggest that the YCS contains different regions that appear to represent transitional regions between the Caribbean and the GoM. According to Olson & Dinerstein (1998), ecoregions are regional-scale conservation units because they encompass similar biological communities, and their boundaries approximately coincide with the area over which key ecological processes most strongly interact. Ecoregional delimitation may affect management actions at regional scales as well as the evaluation of their outcomes (Olson et al., 2001;Edgar et al., 2014). In the case of the large marine system of the GoM, a more precise definition of ecoregions based on the spatial and temporal complexity of the area is required. The inclusion of phytoplankton biomass and species assemblage patterns, in addition to soft-bottom macrofauna, may enhance biological and ecological arguments for defining ecoregional settings in the area. Testing the spatial variation of environmental conditions and biological assemblages along the Southern GoM ecoregion could engender an ecoregional configuration that is more consistent with the ecological variability of marine benthos in the area.