Low Densities of the Ghost Crab Ocypode quadrata Related to Large Scale Human Modification of Sandy Shores

Sandy beaches are the most common ecosystems of coastal regions and provide direct and indirect essential services for millions of people, such as coastal protection, fishing, tourism, and recreational activities. However, the natural habitats of sandy shores are being modified at rates never experienced before, making beaches key monitoring sites of marine ecosystems worldwide. The ghost crab species Ocypode quadrata is the most conspicuous crustacean of sandy beaches along the Western Atlantic coast and has been successfully used as an indicator of anthropogenic disturbance and environmental variability. To investigate the potential role of a “triple whammy” [(1) urbanization; (2) use of resources; (3) decreasing resilience] on the most common bioindicator of sandy shores, we compiled a dataset including 214 records of burrows density from 94 microtidal sandy beach sectors covering a range of over 65° of latitude. The response of burrows density to synergetic effects of human modification of natural systems and environmental changes was investigated using linear models. We used the cumulative Human Modification (HMc) index, a standardized geographic projection of changes of natural systems, as a predictor of urbanization, industrialization and use of resources. The predictor wave energy, tidal range and temperature (sea surface and air) were included as potential effects of climate changes. Literature review showed records mainly concentrated at sub-tropical and temperate regions. HMc values were clearly negatively related to burrows density, thereby supporting an effect of modification of natural habitat at large spatial scale. Sea surface temperature and air temperature were positive related with density and the lack of a general pattern of the relationship between burrows density, interactions between wave energy and tide range, supported unclear patterns reported at regional scales. Finally, we argue that ghost crabs are valuable targets for protection actions on sandy beaches that can benefit coexisting species and provide natural habitat conservation.


INTRODUCTION
Concentrating almost half of the world's population, 1 coastal regions (<100 km from line coast) play a major role in providing ecological services and promote human well-being (Costanza et al., 2014). Although the United Nations 2030 Sustainable Development Goals calls for the protection, restoration, and sustainable use of ecosystems (Cowie et al., 2018), coastal regions are being transformed at rates never experienced before, even in face of evidences that natural systems generate social and economic benefits which exceed those obtained from habitat conversion (Balmford et al., 2002;Xavier et al., 2020).
Recognizing the importance of these services, numerous international collaborations (Benson et al., 2018;Canonico et al., 2019) have centered their efforts to enhance monitoring, setting essential physical, chemical, and marine biological variables to be measured through standardized methods and best practices of data acquisition to subsidize conservation, management, and stackeholders decisions (Pereira et al., 2013;Miloslavich et al., 2018;Fanini et al., 2020;Guerra-Castro et al., 2020). In this sense, sandy beaches are key monitoring sites. Sandy shores are the most common ecosystems of coastal regions and provide direct and indirect essential services for millions of people, such as coastal protection, fishing, and recreational activities (Micallef and Williams, 2009;McLachlan and Defeo, 2018). Negative synergetic effects on biodiversity and ecosystems functions due to multiple human activities co-occurring in a region, are commonly higher when compared to single effects (Raiter et al., 2014;Costa et al., 2020a). This is what Defeo and Elliot (2020) characterized as "triple whammy" on worldwide sandy shores subject to three major threats: (1) increasing urbanization and industrialization; (2) increasing use of resources; (3) increasing susceptibility and decreasing resilience and resistance to the effects of climate change. Thus, a pragmatic implementation of ecosystems management and monitoring efforts on sandy beaches is a major challenge . A complex scenario is presented because sandy beaches should be considered as social-ecological systems, with interactions between human activities and environmental processes involved (Micallef and Williams, 2009;Xavier et al., 2020). Values and conservation goals should be carefully addressed for a broad and successful inclusion of multiple segments of society (Harris et al., 2014;Schlacher et al., 2016).
Monitoring single species, particularly "charismatic" ones due to a public appeal, is perhaps the most common and effective way to monitor environment variability and include social interactions on conservation activities (Micallef and Williams, 2009;Siddig et al., 2016;Xavier et al., 2020). Ghost crabs (Crustacea: Ocypodidae) are the largest crustaceans from sandy beaches with a distribution from the tropics to temperate latitudes (Sakai and Türkay, 2013). Ocypode quadrata (Fabricius, 1787) is the unique ghost crab species from Western Atlantic. Commonly abundant in sandy shores, they are relatively large semi-terrestrial species that inhabit the intertidal-dune interface. Ghost crabs not only respond to anthropogenic drivers, but 1 https://www.un.org/en/about-un/ also to physical environment features such as morphodynamic characteristics (Pombo et al., 2017;Gül and Griffen, 2019) or biological changes like food supply (Tewfik et al., 2016). There are a number of evidences to a common local or regional pattern of burrow density reduction on urban beaches affected by trampling, vehicle traffic, beach cleaning, litter pollution, recreational harvesting, prey depletion and/or armoring (de Souza et al., 2017;Gül and Griffen, 2018;Pombo and Turra, 2019;Costa et al., 2020c). Different preferences of ghost crabs species for across shore zones according to local beach characteristics (Gül and Griffen, 2018;Ocaña et al., 2018) could reveal distinct response to natural habitat modification. Because they are typical upper-shore semi-terrestrial crustaceans, that are theoretically less affected by swash climate and wave action than other macrofauna, they are commonly neglected as a component of beach communities in classical ecological hypotheses (Defeo and McLachlan, 2005;Defeo et al., 2017).
Although some biases of quantification have been raised (Pombo and Turra, 2013), the indirect estimates of ghost-crab population sizes are a non-destructive, friendly, and low cost sampling and, has been used as a valuable and efficient strategy tool by researchers for impact assessments on sandy beaches and today these data represent an important historical of records. At the same time the increasing possibilities to apply remotesensing technologies enable the application of spatially explicit information on conservation and management effort (Monsarrat et al., 2019;El Mahrad et al., 2020). A big effort is being dedicated to make friendly available global scale data set of environmental variables (Fick and Hijmans, 2017;Assis et al., 2018) and intensity of human modification of natural systems, for example the Global Human Footprint (Venter et al., 2016) or the cumulative Human Modification (HMc) (Kennedy et al., 2019(Kennedy et al., , 2020. Associating remote sensing data to the available records of the single ghost crab from Western Atlantic, enabled investigate large scale correlations between environmental drivers and standardized metrics of the magnitude of human modifications of sandy shores. For this purpose, we compiled a dataset including records of burrow density, as a proxy of population density, from pristine to urbanized beach sectors. For the first time, large scale patterns of population densities were related to key marine variables of species distributions and to a standardized geographic projection of changes on natural systems. We brought new insights on O. quadrata distribution and verified the potential role of a "triple whammy" from (1) synergetic effects of resource exploitation and (2) coastal urbanization (HMc), as also from (3) climate changes (wave energy, tidal range, and environment temperature) on the most common bioindicator of sandy shores.

Ghost Crab Data
Studies about the ghost crab O. quadrata were first obtained by searching the databases of Web of Science, Scopus, and Google Scholar using "Ocypode quadrata" as key-word. Additional studies were then sourced by cited references on this pool of documents. For quality control purposes, only studies with peerreview or an equivalent process completed were kept. There is a great variability in methodological approaches used in the literature (for example regarding the area of the beach that is sampled or the metric that counts are reported) that represent a big challenge to investigate spatial and temporal patterns. From this primary source we selected only reports that estimated population size using density of burrows (ind/m 2 ) or all information necessary to estimate this metric. Only studies where densities were estimated from the entire acrossshore area (waterline to upper supralittoral) were used to avoid bias in the estimate of densities. Data that were only available in figures were extracted using the WebPlotDigitizer Program (Rohatgi, 2013 (Figure 1).

Environmental Variables: The Potential Role of Climate Changes
Temperature plays a key role in ghost crab biology because it influences critical physiological and metabolic processes Full, 1994, 1998), survival, distributional limits (Schoeman et al., 2015), reproduction habits (Negreiros-Fransozo et al., 2002) and behavior (Lucrezi and Schlacher, 2014). Mean sea surface temperature (SST, • C) was extracted from the Bio-ORACLEv2.0 database, which offers averages data for the period of 2000-2014. Data source is from Global Observed Ocean Physics Reprocessing (resolution: 0.258/33 vertical levels), and the pre-processed global ocean re-analyses combining satellite and in situ observations generated global data with a resolution of 5 arc-minutes (Tyberghein et al., 2012;Assis et al., 2018). 2 Average values ( • C) of air temperature were obtained from the WorldClim database, which provides long-term bioclimatic data from land areas with a resolution of 5 arc-minutes (Hijmans et al., 2005). 3 The physical environment of sandy beaches is driven by the interaction between wind, sand, waves and tides where resident sandy beach fauna species are thought to be influenced by these physical drivers (Defeo et al., 2017;McLachlan and Defeo, 2018). Data of wave energy (kW/m) were retrieved from the Marine Socio-Environmental Covariates -MSEC database (Yeager et al., 2017). 4 The wave energy values are a product of a wave forecasting model based on wind input for a span of 31 years , including sheltered-coastline corrections (Yeager et al., 2017). The annual average tidal amplitude (cm) was extracted from the AVISO+ data service, which contains global tidal variables with a resolution of 3.75 arc-minutes (Vestbo et al., 2018). 5 The geographical coordinates from each beach sector were used to extract the variables' values. Layers from Bio-ORACLEv2.0, WorldClim and MARSPEC were accessed through the sdmpredictors package (Bosch, 2018) of the R Program (R Core Team, 2020), and layers from AVISO+ were downloaded from the above-cited websites, and manipulated using the raster package (Hijmans et al., 2016).

Human Modification: The Potential Role of Synergetic Effects of Coastal Urbanization and Resource Exploitation
To quantify human modification of natural habitat we used the HMc values (Kennedy et al., 2019(Kennedy et al., , 2020. The index is a degree of human modification across global terrestrial lands and was calculated as the per pixel (1 km 2 ) product (HMs) of the spatial extent and the expected intensity of impact across five major groups of human stressors: (1) human settlement (population density, built up areas), (2) agriculture (cropland, livestock), (3) transportation (major roads, minor roads, two tracks, railroads), (4) mining and energy production (mining, oil wells, wind turbines), and (5) electrical infrastructure (powerlines, nighttime lights). Thus the HMc capture two of the three whammys defined by Defeo and Elliot (2020). For each stressor the median and mean values are from 2016 and 2014, respectively. The final human modification model is defined as: This fuzzy sum is a function that assumes the contribution of a given factor decreases as values from other stressors co-occur. In this way, the HMc is a continuous parsimonious gradient of modification effect that calibrates landscape impacts and ultimately converges to 1.00 (Kennedy et al., 2019).
We also independently accessed the magnitude of human modification of natural habitat at a local spatial scale by applying a categorical classification of impact on the beach record. The classification was primarily assigned following the original descriptions provided by the authors of each study, and calibrated across sites. In some cases, the intensity of environmental modification was not directly declared by the authors; thereby we used Google Earth images and additional public information (i.e., government reports) from each location to assign a common classification. This classification Google Earth images was done based on the image from the year of survey or as close as possible. Our classification included four levels: High Impact -records on beach sectors with completely modified and suppressed supralittoral zone, with easy access (e.g., presence of runway), buildings and markets, and absence of dune/vegetation; Moderate Impact -urbanized beach sectors, but with some degree of conservation of the supralittoral zone and dune vegetation; Low Impact -beaches with difficult or restricted access, relatively preserved supralittoral zone and exclusively impacted by low number of beachgoers; None Impact -totally preserved beaches commonly located at protected areas. Finally, we assessed the relationship between the HMc values and our categorical classification of local human modification. By this way we can estimate some mismatch between the spatial and temporal variability of both indicators.

Data Analysis
For investigate the potential effects of a "triple whammy" on populations of the ghost crab O. quadarta we used linear models.
The density of burrows was included as a response variable. HMc represented the effect of natural modification and the predictor sea surface temperature, air temperature, wave energy and tidal range proxies of the role of climate changes. In a first step a global linear model was run using the response variable burrows density and the linear predictors HMc, SST, air temperature, wave energy, and the interactions wave SST × air temperature and wave energy × tidal range. A fourth-root transformation was applied on burrow densities and the predictors SST, air temperature, wave energy and tide range were standardized. The presence of multicolinearity on predictors was assessed by the variation inflation factor (VIF) and we used a threshold of VIF ≥ 3 to variable exclusion. We applied the Moran's I test (Cliff and Ord, 1981) for spatial autocorrelation in residuals and rejected the null hypothesis of independence. Thus we used the Moran eigenvector approach (Dray et al., 2006) for construct a matrix of spatial patterns that represents the spatial dependence present in the residuals and can be allocated into the linear model. For construct the spatial model we used the package "spdep" (Bivand and Wong, 2018) from the R Program (R Core Team, 2020) and defined a spatial matrix on the indices of points belonging to the set of k = 5, 10, 15, and 20 nearest neighbors. The neighbors of region points were identified by great circle distance and weights matrices were generated. The package spatialreg (Bivand and Piras, 2015) was used to select the optimal subset of eigenvectors and remove spatial dependence on errors. This spatial eigenvectors matrix was included in the linear model and a second stage of model test was applied. Because the model captures the spatial pattern it represents also local environmental variables that were not explicitly inserted in our model. For the validation of the best linear model we plotted residuals against fitted values and each explanatory variable to identify patterns, and violation of homogeneity dependence; histogram of residuals was used to verify normality assumptions (Supplementary Material; Zuur et al., 2009).

Dataset Coverage and Properties
Data included records from Brazil (77%), United States (15%), Mexico (4%), and Cuba (4 distributed from 35.83 • N to 30.28 • S (Figure 1). Sampling year ranged from 1991 to 2019 with a median value at the year 2012 (only one study occurred before the year 2000). From the total number of beach sectors included in the dataset, 81% have temporal replication (at least two records from the same beach sector) ranging from 2 to 5 replicates.
The coastal region analyzed were represent by microtidal sandy beaches of sub-tropical and temperate shores and the summary statistics of the environmental variables are presented on Table 1. HMc values included pristine to very urbanized sandy beaches and varied between 0.036 and 0.92 with a median value equal to 0.53 (Figure 2). We found a general congruence between HMc values and the local impact classification applied. Some degree of mismatch between both indicators revealed that human modification can be captured at distinct spatial scales and possibly temporal and that local characteristics are potential small scale patterns that influence density.

Relationship Between Burrow Density, Human Modification, Environmental and Spatial Predictors
Burrows density ranged from 0 to 1.53 ind/m 2 with a mean of 0.17 ind/m 2 . We identified values of VIF ≥ 3 for SST and air temperature and run two independent models selection including each predictor at each run. Independent spatial matrices were constructed for each model. The best model (using SST) retained the predictors HMc, tidal range, wave energy, SST the interaction tidal range × wave energy and nine spatial eigenvectors ( Table 2), thus supporting the effect of a triple whammy on ghost crabs populations. The best spatial matrix was defined using k = 5 nearest neighbors and we got a significant increase in the model fitting when including the spatial eigenvectors (LR test F = 29.026, p < 0.001). For a model including air temperature this variable was positively related with burrow density (estimate = 0.022) with a marginal significance (p = 0.057) (Supplementary Material). Figure 3 shows the mean partial effects of each predictor on burrows density of O. quadrata. Density clearly decreased with increasing HMc values. We found a positive relationship between burrows density and SST. The interaction between tidal range and wave energy revealed distinct effects along these gradients. We found a positive mean partial effect of wave energy across average and minor values of tide, and a negative effect on beaches with higher values of tidal range.

DISCUSSION
Our results broadly support existent data at local and regional scales indicating a general negative anthropogenic impact on populations of O. quadrata. Using records from sandy beach sectors with multiple morphodynamic states across a range of over 65 • of latitude along the Western Atlantic coast, we showed evidences of a potential role of a "triple whammy" on densities of the ghost crab O. quadrata. Herein density was clearly negatively related to HMc values that capture the extent of human settlement measured by metrics such as population density, builtup areas and the presence of major roads. For the first time a global standardized metric of human modification of natural system was related to densities of a sandy beach species, thereby supporting the effect of large scale urbanization (sensu Defeo and Elliot, 2020) on the biodiversity of sandy shores. Environmental variables were also important predictors of burrow densities and, the interaction between tidal range and wave energy effects evidenced the lack of a general relationship between ghost crab distributions and environmental drivers. Local spatial variability could only be captured by the spatial model due to our spatial resolution of the environmental predictors. This potentially include unknown precitors explaining density at those scales. Because the same or neighbour sectors were sampled at different times this spatial model could be, in some way, biased by a temporal effect.
The main human intervention on sandy shores consists of physical modifications which have critical effects on coastal ecosystem dynamics by reducing natural habitat suitability (Jones et al., 2017). Particularly higher human population sizes are related to easy access, efficient public transportation and recreational facilities as restaurants, bars and hotels over beach sectors that enhance their social value. In fact, this is the case of hundreds of beaches from the coast of Brazil (Amaral et al., 2016). Recreation-associated disturbances, such as trampling and vehicle traffic, clearly impair invertebrate macrofauna (Gheskiere et al., 2005;Schlacher and Thompson, 2012), including ghost crabs (Schlacher et al., 2007;Lucrezi et al., 2009;Costa et al., 2019Costa et al., , 2020c. The suppression of the dry upper beach zone extirpates semi-terrestrial invertebrates due to habitat loss Hubbard et al., 2014;Cardoso et al., 2016). Although ghost crabs can cover the entire extent of the beach face, they depend directly on supralittoral to mid-intertidal zones to feed and to construct their semi-permanent burrows (Lucrezi and Schlacher, 2014).  The levels in the categorization used here to measure local impact relied mainly on physical modification of natural features characterizing the urbanization effects of the "triple whammy." Some mismatch between this classification and HMc values revealed important distinctness between the spatial scales of the metric to be used. For example, we can have a preserved supralitoral zone of a beach sector located at a metropolitan area. Backshore vegetation is successfully used as a metric of natural habitat modification of sandy beaches (Schlacher et al., 2008;Orlando et al., 2020) and is an important habitat for O. quadrata, providing shelter, burrow stabilization, foraging area, and protection from storms (Lucrezi et al., 2010;Schlacher et al., 2011;Gül and Griffen, 2019). This is also true for other species of macrofauna. Cardoso et al. (2016) reported that despite having some infrastructure facilities, beaches with well-preserved backshore vegetation in Rio de Janeiro state, Brazil, supported higher densities of the sandhopper Atlantorchestoidea brasiliensis, a mobile crustacean typically found from the supralittoral to the upper midlittoral zone (Cardoso and Veloso, 2001;Gómez et al., 2013). Similarly, Orlando et al. (2020) reported that macrofauna community richness is usually higher in areas with preserved vegetation on sandy shores. Therefore, it is expected that ghost crab conservation confers a "protection wave" to numerous coexisting species (beneficiary species "under the umbrella") (Simberloff, 1998;Fleishman et al., 2001), and can co-occur with the control of erosive processes through the maintenance of coastal vegetation.
Although neglected in sandy beach impact assessments (González et al., 2014;Luarte et al., 2016) nightlight, also weighted on HMc values and related to urbanization, can affect ghost crab sensorial responses (Silva et al., 2017). Costa et al. (2020b) tested whether artificial light attracts ghost crabs inland toward roads nearby sandy beaches, but evidences that light pollution predicts ghost crabs road-kills (car-crab collisions) were not found. In an extensive meta-analysis, Schlacher et al. (2016) pointed out that mechanistic processes and suitable sampling design concerning impacts of artificial nightlight on macrofauna species remain underreported, even though foraging behavior and activity rhythms may be sensitive to artificial light and could possibly have fitness implications for ghost crabs species.
The resilience to synergetic effects of these rate of coastal modification and environmental change includes the combination of the capacity to resist to increasingly frequent and severe disturbances and to adapt to new environmental conditions, that is enhanced by high population sizes and dispersal rates (Bernhardt and Leslie, 2011). Thus, we expect that population decrease on urbanized sandy shores will increase susceptibility and decrease resilience and resistance to the effects of environmental changes. Understand how these changes will disturb ghost crab populations includes establish ecological theories about patterns of macrofauna communities of sandy beaches to morphodynamic features.

The Swash Exclusion Hypothesis and the Habitat Harshness
Hypothesis propose that macrofauna along a morphodynamic gradient is limited by the swash climate; and then, species richness, abundance and biomass increase from reflective to dissipative beaches (Mclachlan et al., 1993;Defeo and McLachlan, 2005). For our data the interaction between tidal range and wave energy indicated the importance of these environment predictors to population densities. A positive effect of wave energy was estimated to sandy beaches with lower values of tide range, particularly characteristic from temperate regions of South America with dissipative beaches with high wave energy and high primary productivity (Odebrecht et al., 2014). A weak but negative effect of wave energy occurred for higher values of tide range, that are characteristic of sandy beaches from southeast coast of Brazil (Amaral et al., 2016), the region that included the majority of records in our data. In ghost crabs literature we have quite different results. For example, Quijón et al. (2001) reported no clear response of Ocypode gaudichaudii to morphodynamic variables on sandy beaches in the North of Chile. Counting indistinguishable burrows of three ghost crab species, Lucrezi (2015) observed an increase in burrow densities toward dissipative conditions on four warmtemperate microtidal sandy beaches in South Africa. Oppositely, Defeo and McLachlan (2011) found increase, but not very clear, in burrow densities of O. quadrata from dissipative to reflective beaches and coarser grain sizes, compiling information from microtidal sandy beaches from the warm temperate Southwestern Atlantic province. Finally, analyzing densities along nine pristine areas with distinct morphodynamic states and wave exposure levels, Pombo et al. (2017) reported no clear response to slope, but an increased burrow density in finer-grained beaches. These evidences supported that beach morphodynamic state alone could be insufficient to predict population and community parameters, that impose some challenges in proposing forecasting scenarios, especially for supralittoral species which are not directly exposed to swash conditions (Contreras et al., 2003;Defeo and Gómez, 2005).
Because O. quadrata has an indirect development with a long-lived larval stage and populations sharing haplotypes over ∼7,000 km (Mattos et al., 2019) is expected that not only coastal modification and beach characteristic will influence population densities, but also oceanic variables. Although sandy beaches occupy about one-third of ice-free coastline (Luijendijk et al., 2018), large scale (i.e., regional and continental) oceanic environmental drivers of macrofauna distribution have only been investigated quite recently (Defeo et al., 2017). The positive effect of SST, as well as of for air temperature when SST is excluded, should be evaluated with care because it seems that these predictors can act distinctly on ghost crabs distribution according to ontogenetic stage. O. quadrata has a planktonic larvae phase, thus sea temperature is more likely to be related to larvae metabolism, development rates, metamorphosis and time spent in the plankton (O'Connor et al., 2009;Gerber et al., 2014). Recently, Meerhoff et al. (2020) showed the influence of seasonal winds and associated ocean circulation patterns to larval connectivity of the mole crab Emerita brasiliensis. Local hydrodynamic regimes may have great potential to influence larval settlement. On the other hand, air temperature is likely to be related with population establishment, distribution range and burrowing activity. Schoeman et al. (2015) showed that distribution range of Ocypode cordimanus is not primarily determined by recruitment limitation along the Central to Northern East coast of Australia. In that region, both air and ocean temperatures decrease from north to south, but with a stronger spatial gradient in seawater than in air temperatures. Abundance of species is usually lower near its range edge (Bates et al., 2014) and this is in part coherent with our findings, which show decreasing burrow densities toward lower sea surface and air temperatures.
In conclusion, we showed evidences that coastal modification plays a potential major role on the population density of the ghost crab O. quadrata along the Western Atlantic coast and, these effects should vary across environment characteristics. Effects at continental scales seems not only disturb densities of but also individuals size (Schlacher et al., 2016). The species is a model of key element in sandy beach trophic webs and thus indicates that natural habitat loss at a continental scale will certainly bring changes on patterns of the land-sea interface processes. Urbanization of sandy beaches has an intrinsic relation with widely used indicators of natural habitat urbanization, such as human population density and nighttime light intensity, here identified as potential and promising proxies of biological disturbance at large spatial scales. Ghost crabs conservation will benefit coexisting species, because it requires protection to sandy shores' natural habitat. In addition, maintenance of coastal vegetation was found to be essential for ghost crabs conservation, and it indirectly allows the control of erosive processes. Therefore, targeting ghost crabs as umbrella and indicator species is a valuable strategy for the conservation and monitoring of coastal regions. Our dataset is may be influenced by intrinsic processes of the Brazilian coast, the area with the greatest research effort for this species. Also, data is restricted to microtidal sandy beaches with a lack of existent reports from tropical coasts and macrotidal regimes. However, we argue that the number of records included here, with a wide range of morphodynamic states across a continental scale covering the entire gradient of urbanization levels, supports our findings in this domain and, brings new insights for the ecology of ghost crabs. We strongly recommend increasing monitoring data on the tropical coasts of the Western Atlantic, mainly sandy beaches along the Caribbean Sea and the Gulf of Mexico. Finally, manipulative experiments using standardized protocols that include key environmental variables at large scales are pivotal, particularly to disentangle causality relationships between reduction of ghost crab populations and human stressors.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CB led the analyses and writing of the manuscript in close collaboration with GM and LC. AS-G and IZ contributed critically to the drafts and gave final approval for publication. All authors contributed to the article and approved the submitted version.  (CNPq Proc. 301203/2019-9). AS-G was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq Proc. 301475/2017-2). LC was supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES Proc. 88882.463168/2019-01).