Nestedness patterns and the role of morphodynamics and spatial distance on sandy beach fauna: ecological hypotheses and conservation strategies

Sandy beach fauna is hypothesized to be mainly structured by environmental variables. As such, it is expected that morphodynamic characteristics are limiting factors, and the species pool inhabiting harsher reflective beaches would be a subset of (i.e., nested in) the fauna of nearby dissipative beaches. We investigated the existence of a nestedness pattern in sandy beach assemblages, as well as the contribution of environmental and spatial variables (i.e., factors that potentially affect an assemblage regardless of environmental conditions - typically related to distance between sites and dispersal of organisms) on sandy beach macrobenthic fauna. Dissipative beaches had higher species richness than reflective beaches but we found no nestedness pattern. Furthermore, almost every beach showed exclusive species. Spatial variables exerted stronger influence on macrobenthic assemblages than local environmental variables. Our results therefore suggest that local and small-scale recruitment is the predominant process structuring macrobenthic assemblages. These results bring important implications for sandy beach conservation: given that spatial distance is an important factor structuring macrobenthic fauna and different sandy beaches harbor different pools of species, conservation programs need to focus on sandy beaches across large spatial scales and with varied morphodynamic characteristics in order to preserve coastal biodiversity.

(1) Given that sandy beach fauna is tightly related to morphodynamic characteristics, environmental variables should be the main structuring factor of intertidal macrobenthic assemblages. Furthermore, at a regional scale, beaches with similar morphodynamic characteristics should harbor similar macrobenthic assemblages. (2) Considering that the harsher environmental conditions towards reflective beaches determine the progressive exclusion of species, the species pool of a reflective beach would be a subset of species from nearby dissipative beaches. In other words, at a regional scale, biological communities from reflective beaches would be nested in the species pool of dissipative beaches. (3) Finally, for a nestedness pattern to arise from beaches with different morphodynamic characteristics, it is necessary that dissipative beaches not only have higher species richness, but also a higher number of exclusive species than reflective or intermediary beaches.
Recently, a few studies have been investigating some of those assumptions, but a comprehensive evaluation of the main hypotheses on distribution patterns of sandy beach macrobenthic communities is yet to be conducted. The influence of environmental and spatial variables (i.e., factors that affect a local assemblage regardless of environmental conditions -typically related to distance between sampling sites and dispersal of species) on sandy beach assemblages has been examined 19,20 ; however, to our knowledge, only one attempt was made to estimate nestedness patterns in sandy beaches 21 . Since the publication of this study, a significant advancement has been made in statistical procedures to estimate nestedness, which allows for more robust calculations of distribution patterns 22 .
Knowledge about nestedness patterns on sandy beaches is important not only to improve our understanding of their ecological functioning, but also to provide valuable information for biodiversity conservation and management strategies. Nestedness patterns inform us on how species are distributed across different habitats 23,24 , information that is crucial to identify areas that should be prioritized for protection 25,26 . A nested pattern in sandy beaches would imply that the set of species inhabiting reflective beaches is a subset of the set of species inhabiting dissipative beaches; therefore, dissipative beaches should be priority targets for biodiversity conservation. On the other hand, a non-nested pattern, especially with exclusive species occurring in dissipative and reflective beaches, would suggest that a combination of different types of sandy beaches should be preserved. This information is crucial for sandy beach ecosystems, which are often neglected in coastal conservation programs 16,27 .
The aim of this study is twofold: (1) to improve the understanding of macrobenthic distribution patterns on a regional scale by comprehensively testing sandy beach ecological hypotheses, and (2) to provide information that can be used in conservation and management of coastal habitats. Specifically, we investigated (i) the influence of environmental and spatial variables on the distribution of sandy beach macrobenthic assemblages, (ii) the dissimilarity of macrobenthic assemblages in a set of nearby sandy beaches with contrasting morphodynamic characteristics; (iii) their nestedness patterns and the occurrence of exclusive species.

Results
Environmental characterization and biological data. Sampled beaches/sectors showed a range of morphodynamic characteristics (Table 1). According to Beach Index (BI) values and beach width, Toque-Toque and Picinguaba were characterized as reflective; Barra do Sahy, Baleia, Flecheiras 1 and 2, Cidade 2 to 4 and Camaroeiro were characterized as intermediary; whereas Cidade 1, Fazenda 1 and 2 and Palmeiras showed dissipative features. There was very little variability in BI values between periods. Beach slope and mean diameter followed the trend shown by the beach index, with reflective beaches having coarser and dissipative beaches having finer grains. Beach width varied greatly, ranging from 25 m (Picinguaba) to 140 m (Baleia and Fazenda 2). Calcium carbonate content also showed a great variation, ranging from 0.35% (Picinguba) to 7.33% (Baleia). Organic matter content (%) was generally low across the studied beaches (0.13% at Sahy to 1.39% at Palmeiras beach).
A total of 2834 individuals, comprising 112 macrobenthic species, were collected throughout the study. Mollusks contributed to 53.6% of the total number of species (42 species of gastropods and 18 species of bivalves), whereas polychaetes accounted for 35.7% (40 species), crustaceans for 8% (9 species), and sipunculans and ophiuroids accounted for 1.8% each (2 species each). Mollusks and polychaetes also accounted for more than 95% of macrobenthic abundance, with a slightly higher number of mollusks being sampled (1440 and 1297 individuals, respectively).
Influence of spatial and environmental variables on sandy beach macrobenthic fauna. Beach width, slope and mean diameter were the best explanatory environmental variables of macrobenthic distribution. Selected eigenvectors (dbMEM) were mainly related to large-scale patterns ( Community dissimilarity among beaches. Dissimilarity in species composition of macrobenthic communities seemed to be more related to spatial distance than beach type. Macrobenthic assemblages were more similar among intermediate beaches, with the exception of Sahy and Baleia, the two intermediate beaches located further away from the others. Fazenda beach sectors (characterized as dissipative) had a stronger similarity with the neighboring beach, Picinguaba (reflective), than with the dissipative beach Palmeiras. Beach type, however, was an important factor, as Palmeiras beach was different from the neighboring intermediate beaches, and Toque-Toque, the only reflective beach on the southern part of the study area, showed a high dissimilarity with other beaches (Fig. 1).  Nestedness patterns. NODF (Nestedness metric based on overlap and decreasing fill) values show a non-nested pattern of species distribution along the investigated environmental gradients of each variable. In fact, the observed values were consistently lower than the expected values, with negative z-scores, indicating a tendency towards an anti-nestedness pattern. The differences between NODF observed and expected values, however, were only marginally significant for all ordering variables in both periods (p < 0.10). The results for Matrix Temperature (T) were not significant, and a non-nested pattern was recorded for most periods and variables. With the composition from both periods pooled, the tendency towards anti-nestedness was much lower, and the observed nestedness was not different from the expected value under the null model for every variable and for both NODF and Matrix Temperature (Table 4).
Relationship between number of species and beach type. The number of species was positively correlated to BI (n = 28, r = 0.40, p = 0.037), negatively correlated to beach slope (n = 28, r = −0.46, p = 0.014), but not to grain size (n = 28, r = 0.30, p = 0.126) and intertidal width (n = 28, r = 0.11, p = 0.580). No increase   of exclusive species towards dissipative beaches was observed. In fact, most beaches had 1 to 3 exclusive species registered, regardless of BI value (Fig. 2). Contrary to our predictions, the highest number of exclusive species was registered on intermediate beaches.

Discussion
Significant changes in macrobenthic assemblages were associated with different sandy beach types, but morphodynamic characteristics alone did not explain the observed biological patterns. A higher number of species was found in dissipative beaches, and environmental characteristics significantly explained the structure of sandy beach assemblages, partially supporting sandy beach hypotheses. However, distance between beaches had an overall higher contribution to the variation in species distribution than environmental variables. Furthermore, we did not find a nested pattern in the morphodynamic continuum nor an increase of exclusive species towards dissipative beaches. The number of macrobenthic species was significantly correlated with BI and beach slope, highlighting that environmental characteristics are an important factor structuring intertidal macrobenthic assemblages 4,7,28 . This importance was further reinforced by the significant contribution of the environmental variables alone in explaining the distribution of macrobenthic assemblages. Nevertheless, our results also revealed that the distribution of macrobenthic organisms in the study area was mainly related to spatial variables. Similar results, i.e. strong influence of spatial variables on marine soft sediment assemblages, have been reported in recent studies 19,20,29,30 and are an indication that spatial processes are an important structuring factor in this system.
Several authors have pointed out that the variation accounted by spatial processes may arise from two main sources: it can be attributed to pure spatial processes related to dispersion of organisms or to some spatially structured unmeasured environmental variables [31][32][33] . Here, it is likely that both sources affected sandy beach macrobenthic assemblages.
Sandy beaches are thought to be systems highly connected by species dispersal, forming metapopulations that are limited by environmental conditions 5,16 . Although this process holds true at the population level for some macrobenthic taxa 15,34 , it may not occur at a community level. The stronger influence of the spatial component alone in comparison to the environmental component, and the high similarity of macrobenthic assemblages in neighbouring beaches, despite their morphodynamic state, suggest that local recruitment and small-scale dispersion are predominant processes for macrobenthic metacommunity dynamics at the study area. In this sense, our results agree with suggestions that sandy beaches are somewhat isolated systems 35 (although the degree of isolation is likely to vary among species and dispersal strategies 20 ). Our results also agree with evidence that macrobenthic species tend to recruit locally, despite dispersal potential, due to the favourable conditions of the parental area 36,37 . It is important to notice that intertidal species are more limited in dispersion capability than subtidal ones 38 , and studies evaluating the subtidal community may find spatial processes to contribute differently.
On the other hand, spatially structured unmeasured environmental variables could also be responsible for the stronger influence of the spatial component. Although we included many environmental variables that are commonly found to affect the species distribution of soft sediment macrobenthic assemblages (e.g., slope, width, BI and sediment), many others which could be spatially structured were not considered. Physical forces such as local currents, winds, and small-scale turbulent processes in the water column are known to affect recruitment, distribution and mortality of marine invertebrates [39][40][41][42][43] , and, therefore, may contribute significantly to their spatial patterns. Oxygen and nitrogen content in the sediment, salinity and algal biomass are other important environmental factors that were not considered here 17,19,44 . Thus, it is likely that the environmental control observed here is underestimated and it is important that future studies on macrobenthic assemblages include more variables than only morphodynamic characteristics.
Although the inclusion of a larger set of environmental variables may increase the importance of the environmental component explaining the patterns of richness and distribution of sandy beach macrobenthic assemblages, our results clearly show that spatial processes, which are often neglected in sandy beach ecological studies 4,7 , can exert a significant influence on these assemblages. Here we considered distance as the spatial component; nonetheless, the use of connectivity measures in the spatial structure has been shown as a promising tool in the evaluation of spatial processes in metacommunities 45,46 and should be tested in further studies on sandy beaches.
Nestedness usually arises as a result of a species filtering process along an environmental gradient due to differences in the tolerance of species to that gradient 22,47 . Sandy beaches show a clear environmental variation along the morphodynamic continuum; however, our expectations of significant nestedness patterns on sandy beach macrobenthic fauna were not corroborated, which shows that macrobenthic assemblages of reflective beaches are not a subset of those from dissipative beaches. The higher contribution of spatial processes compared to environmental variables to the distribution of macrobenthic assemblages in the study area is likely to reduce the influence of the species filtering process, resulting in the non-nested pattern observed.
In some cases, the macrobenthic metacommunity presented a tendency towards anti-nestedness. This concept is used in different situations, and thus the origin of such pattern is hard to distinguish 48,49 . One explanation may lay in low prevalence, where species presented in a single site increase the tendency towards anti-nestedness 50 . Exclusive species were found on almost every beach and during every period in the study area. This pattern resulted in a sparsely filled species by site matrix (i.e. with many non-occurrences). When more than 40% of sites within a dataset are required to include all species, these datasets show a tendency towards anti-nestedness patterns 49 . These factors may be responsible for the observed tendency towards anti-nestedness patterns in the distribution of macrobenthic species. In this sense, anti-nestedness can be seen as a process of species replacement and high beta diversity 51 .
Besides increasing the ecological understanding of sandy beach assemblages, our results bring important implications for conservation programs and management plans of sandy beaches. When nestedness occurs, species-rich assemblages can be prioritized for conservation, as species-poor assemblages comprised only a subset of species 52 . Our results, however, show that different beach types are likely to contain complementary sets of species, as showed by the non-nested patterns and the dissimilarity among sites. This complementarity is further reinforced by the presence of exclusive species in almost every beach analyzed. Therefore, sandy beaches with different morphodynamic characteristics need to be protected in order to maintain overall biodiversity, and management programs with this goal should not be based only on sandy beach type.
Overall, we found partial support for hypotheses explaining the distribution patterns of intertidal macrofaunal assemblages. Morphodynamic characteristics were an important factor structuring the fauna, but spatial distance had a stronger influence on macrobenthic distribution. We found no support for the hypothesis of a nestedness pattern related to morphodynamic beach type. Instead, we detected a non-nested distribution, with exclusive species occurring on almost every beach. These results indicate that processes related to dispersal may be more important than local environmental features and that different sandy beaches harbor different assemblage composition. Therefore, we suggest that conservation practices across regional/large spatial scales should consider beaches with distinct morphodynamic features in order to preserve regional biodiversity. Future studies should investigate a more complete set of parameters attempting to find conservation targets that adequately protect sandy beach ecological diversity. This is an imperative task on which our quality of life in the long term relies 16 .

Methods
Study area. This study was conducted at the Northern coast of the state of São Paulo, Southeast Brazil, in the municipalities of Ubatuba, São Sebastião and Caraguatatuba (Fig. 3). The area has a Cfa Köppen's climate characterized as humid subtropical with hot summers and a lack of dry season 53 . Nine beaches in this area were selected for sampling (Barra do Sahy, Baleia, Palmeiras, Flecheiras, Cidade, Camaroeiro, Fazenda, Toque-Toque and Picinguaba). Beaches with a length greater than 1 Km (Cidade, Flecheiras and Fazenda) were divided in sectors due to different features registered along the shore, such as beach width and grain size, that could directly influence macrofauna abundance and occurrence 54,55 . Taking this into account, Cidade beach was divided in four sectors (Cidade 1 to 4), Flecheiras and Fazenda beaches in two (Flecheiras 1 and 2, Fazenda 1 and 2). All other beaches consist of 1 sector only, leading to a total of 14 sectors sampled. The position of each sector was recorded with a GPS (Garmin eTrex Legend, datum WGS84) and the same locations were sampled during each sampling event. Each beach sector was sampled twice, once during austral Fall/2001 and once during Spring/2001.

Sampling desing.
Sampling was done at three sampling sites (area of 100 m² each, 10 m on each side), equally spaced in the intertidal zone from the water line to the drift line, in the central area of each sector. Given the wide range of sizes of macrofaunal taxa, from a few millimeters to tens of centimeters, two complementary methods were used in order to obtain the most accurate estimate of biodiversity. Thus, in each sampling site, five samples were randomly collected with small corers (d = 10 cm, A = 0.008 m²) and sieved on mesh size of 0.5 mm to evaluate small-bodied fauna; and three samples were randomly collected with a large corer (d = 45 cm, A = 0.16 m²) and sieved on mesh size of 1.0 mm to evaluate larger species.
Environmental characterization. To characterize the different sectors sampled, three sediment samples (A = 0.008 m²) were taken at each sampling site to evaluate granulometrical parameters (grain size and sorting coefficient), organic matter and calcium carbonate content (CaCO 3 ). Beach type (dissipative, reflective or intermediate) was classified using the Beach Index (BI) 7 , according to equation (1): 10 where Mz is the mean grain size (φ + 1); TR is the maximum spring tide (m) and S is the beach slope (1/tanβ). Higher values of BI denote more dissipative characteristics.
Influence of spatial and environment variables on sandy fauna. To investigate the relative contribution of spatial distance and environmental variables on sandy beach macrofaunal assemblages, we used a variation partitioning procedure 56,57 applied to the redundancy analysis (partial RDA) 31 . Abundance data from both periods were pooled to increase the number and abundance of species per site. Pooling was justified by Procrustes analysis results 58 , which showed that biological datasets from both periods were strong and significant correlated (Correlation = 0.823, P = 0.037). Sediment mean diameter was also not different between periods at each beach (Supplementary Table 1), given that beach width and slope were constant between periods, beach index (BI) was also not different. Other variables were not used in the redundancy analysis. We used slope, beach width, beach index and mean grain size of each location as environmental variables. To derive spatial variables (proxy for spatial processes), we used distance-based Moran eigenvector maps (dbMEMs) for the matrix composed by the geographical coordinates of the sites. dbMEMs provide orthogonal vectors that maximize the spatial autocorrelation 57 . Large-scale spatial correlation is modelled by the initial dbMEMs, while the last dbMEMs correspond to fine-scale spatial correlation 57 . We used the longest distance connecting two neighboring sites as a threshold to truncate the distance matrix. Distances larger than the threshold value were replaced by an arbitrarily large value equal to four times the threshold and were disconnected in a neighbor matrix (i.e., truncated matrix) 59 . To minimize random effects by dominant taxa and make data more appropriate to be analyzed by linear ordination methods, community abundance data was transformed using Hellinger function 60 The variation was partitioned in four fractions: (a) influence of environmental variables alone (slope, beach width, beach index and mean grain size); (b) shared influence of environment and spatial variables; (c) influence of spatial variables alone (distance-based Moran eigenvector maps, dbMEM); and (d) unexplained variation. The significance of each component was assessed by permutation tests (n = 10000) on the redundancy analysis model. Variable selection for both components was done using forward procedure with double stopping criteria, which considers Akaike values and Adjusted R² to select variables, reducing type I error likelihood 61 . Only variables selected by this procedure were included in the variation partitioning analysis.
Community dissimilarity among beaches. We used a non-metric dimensional scaling ordination (NMDS) to assess the similarity in community composition among sectors. Abundance data for the full suite of species was log-transformed (log x + 1) to down weigh the influence of extreme values. Bray-Curtis index was used to calculate the community composition dissimilarity among areas. Similar to the redundancy analysis, data was also pooled for this analysis. The similarity in community composition and variation partitioning were carried out using the R Software v. 3.3.1 62 , with the additional packages vegan 63 and spacemakeR 64 .
Nestedness analysis. The community matrix of each sector was arranged with presence/absence records for each species. In order to evaluate nestedness, the species-by-site matrix needs to be ranked. The common procedure is to rank species in decreasing order of frequency and sites in decreasing order of richness, but ranking can be made according to a gradient of interest, in order to check whether nestedness arises from changes in the environment 22 . Thus, to test our hypothesis that reflective beaches have a subset of species from dissipative beaches, we ranked sites according to four proxies of beach morphodyamics: (1) decreasing value of Beach Index (BI); (2) decreasing value of mean grain size (φ); (3) increasing value of beach slope (1/S); (4) decreasing beach width (m). This procedure was carried out for individual as well as for pooled samples, in order to account for possible temporal variations in nestedness.
Two different metrics were calculated to estimate nestedness: the matrix temperature (T) 65 and the nestedness overlap and decreasing fill (NODF) 23 . Each metric has a different theoretical basis to estimate nestedness and may show contrasting behavior 23,66 . T is a metric of matrix disorder, accounting for unexpected presences and absences from a perfectly nested matrix, ranging from 0° (perfectly nested) to 100° (minimum nestedness) 65 . NODF is based on the overlapped presence and the marginal totals between pairs of rows and columns, ranging from 0 (no nestedness) to 100 (maximum nestedness) 23 . The use of more than one metric to quantify nestedness is a common practice 49,67 , and we used the two metrics to check if an observed pattern is consistent among metrics. Simulations were made under the FF null model (Fixed-Fixed, n = 1000 randomizations), which randomizes preserving row and column total value. This null model has a more conservative approach (i.e., less prone to type I error) and also has a low sensitivity to matrix fill, size and shape 22,23 . Z-scores were calculated as a mean of standardization of values, following equation (2): where X is the observed index value, µ is the mean and σ is the standard deviation of simulated expected index values. For NODF, positive values indicate a tendency towards nestedness, whereas negative values indicate a tendency towards anti-nestedness. For T, the opposite trend is expected. Nestedness metrics and null model analysis were carried out using the NODF software 68 .
Relationship between number of species and beach type. The relationship between the number of species and the number of exclusive species with the Beach Index was carried out using Pearson linear correlation. Periods were pooled for this analysis to give a more accurate estimate of macrobenthic richness.
Data availability. The datasets generated and analyzed during the current study are available upon request to the corresponding author.