The study area comprised 12 sandy beaches along the Spanish coast of the Gulf of Cadiz, from Hoyo beach (37°11’55”N, 07°17’45”W) near the Portuguese border to Lances beach (36°02’31”N, 05°38’08.031”W) in the area near the Strait of Gibraltar (Fig. 1).

For all 12 beaches, macroinfauna sampling was carried out in six replicated transects haphazardly separated by a 100-m-long stretch in the central area of each beach during low tide between March and May 2011.

In each transect, macroinfauna samples were collected at 10 equally spaced shore levels extending from above the drift line (level 10) to the low swash zone (level 1) with a 25-cm-diameter plastic core of a depth of 20 cm. Transect length depended on each beach’s width. A total of 60 samples were collected within a total sampled area of 3.75 m² per beach. On temperate beaches this area is considered sufficient to collect 90% of all the macroinfauna (Jaramillo et al. 1995). Samples were sieved on site through a 1-mm mesh sieve, collected in labelled plastic bags and preserved in 70% ethanol stained with Rose Bengal. Additionally, one sediment sample was taken at each sampling level with a plastic tube (3.5 cm diameter) buried 15 cm deep to analyse the mean grain size, sorting coefficient (Trask 1950), sand moisture and organic matter of the sediment.

In the laboratory, the macroinfauna was quantified and identified to the lowest taxonomic level possible. The mean grain size was determined following the method proposed by Guitián and Carballas (1976). This method discriminates different granulometric fractions when the sediment composition is mainly sand and the pelitic fraction is low (less than 5%). Sand moisture was determined by measuring the weight loss after drying the samples at 90°C. The organic matter content of the sediment fraction <0.5 mm was estimated as the difference between dry sediment weight and sediment weight after calcination at 500°C.

To characterize the morphodynamic state, the relative tidal range (Masselink and Short 1993), the beach index (McLachlan and Dorvlo 2005), the beach state index (McLachlan et al. 1993) and the dimensionless fall-velocity parameter (Dean’s parameter; Dean 1973) were used. The beach face slope was estimated by the height difference according to Emery (1961). The height and wave period, temperature and tidal range were taken from average annual data from the SIMAR oceanographic database of Puertos del Estado (Spanish Ministry of Public Works) (http://www.puertos.es/es-es/oceanografia/Paginas/portus.aspx).

To characterize the main use of beaches, the conservation index (CI) and the recreation potential index (RI) were applied following McLachlan et al. (2013). For each indicator a simple scoring system was developed based on three key drivers in order to evaluate the potential/suitability of a beach for conservation and/or recreation. CI takes into account 1) the extent, nature and condition of the dunes, their well-developed vegetation and their connection with the beach; 2) the presence of iconic and endangered species; and 3) the macrobenthic community abundance and species richness. In contrast, RI is based on 1) available infrastructures to support recreational activities (beach access, toilets, etc.); 2) beach safety and health status; and 3) physical carrying capacity. CI and RI values range from 0 to 10, in increasing order of conservation value or recreation potential.

Univariate analyses were used to characterize the fauna communities present on each beach studied, calculating the Margalef richness index (d), the Shannon-Wiener diversity index (H’) and the Pielou evenness index (J’), using the PRIMER software package.

The zonation pattern on each beach studied was identified using cluster analysis, based on the Bray-Curtis similarity matrix, followed by a similarity profile test (SIMPROF) (Clarke and Gorley 2006), to evaluate the significance of the classification (p<0.05). Previously, abundance data were fourth-root transformed to downweight the contribution of the major abundant species.

Once the zonation patterns had been defined on individual beaches, a modal pattern of zonation was established for the entire set of beaches. To this end, species from each sampling level were pooled based on zones identified by cluster analysis. Then, a single matrix of “species × zone” for each beach was generated and all of them were combined into a global matrix. This global biological matrix was fourth-root transformed and subjected to non-metric multi-dimensional scaling ordination (n-MDS). Furthermore, the similarity percentages analysis (SIMPER) was performed in order to find the typifying species in each identified zone for the entire set of beaches along the Gulf of Cadiz. Beaches that did not have a clear zonation pattern were excluded from these analyses. All multivariate analyses were performed with PRIMER-E v 6.1 (PRIMER-E Ltd) (Clarke and Warwick 2001).

To determine associations of macroinfauna communities with environmental variables, a canonical correspondence analysis (CCA) was applied. First, a global biological matrix was submitted to detrended correspondence analysis (DCA) in order to measure the gradient lengths and to ensure an unimodal species response. Gradient length of the first axis was greater than 3.0 SD, and a CCA ordination method was used. For this analysis, only the most abundant species (>6% of total contribution in each biological zone identified) were taken into account after fourth-root transformation.

The environmental parameters matrix was transformed (Log(x+1)) and standardized prior to reducing extreme values and providing better canonical coefficient comparisons. Only variables significantly related to the fauna variation were included (p<0.05); to this end, each variable was analysed separately and its significance was tested using a Monte Carlo permutation test (999 permutations) (TerBraak 1995).

In CCA analysis, the statistical significance of canonical eigenvalues and the significance of the first two axes were tested by a Monte Carlo test (999 permutations). DCA and CCA were carried out with the statistical software package PC-ORD (McCune and Medford 1997).

The physical and morphodynamic characteristics of the 12 beaches studied are shown in Table 1. According to the tidal range and relative tidal range, the beaches were categorized as mesotidal, dominated by waves. The beaches showed a wide range of morphodynamic types, classified by Dean’s parameter as intermediate (Terrón, Mazagón, Matalascañas, Barrosa and Zahara), dissipative (Bota, Costa Ballena, Valdelagrana, Levante, Cortadura, and Lances) and reflective (Hoyo). Beach state index values classified most beaches as intermediate to dissipative, with high energy levels, except for Hoyo and Zahara, which were intermediate beaches with low to medium energy levels.

A total of 63 macroinfauna taxa were recorded from the beaches of the Gulf of Cadiz (Appendix 1). Crustaceans were the most diverse taxa, with 23 species, followed by polychaetes (22 species), insects and molluscs (9 and 8 species, respectively). Appendix 1 shows the total abundance, total species, Margalef species richness, Shannon-Wiener diversity index and Pielou evenness index. Bota and Levante had the highest richness, with 33 and 24 species, respectively, while the lowest value was found at Matalascañas (4 species). Density was also highly variable, ranging from 23 to 249 individuals per m². The lowest diversity value (H’) was observed at Matalascañas beach (0.40), while the highest was found at Levante (2.68). The evenness index ranged from 0.29 at Matalascañas to 0.86 at Lances beach.

Scolelepis spp., Pontocrates arenarius and Haustorius arenarius were the most frequent species. The polychaete was present on 100% of the total beaches sampled, while the amphipods were present on 90%. In terms of abundance, polychaetes of the genus Scolelepis were dominant, accounting for 28% of total abundance, followed by the amphipod Haustorius arenarius, which accounted for 15% of the total.

Across-shore species distribution on each beach studied is shown in Figure 2. Cluster ordination and the SIMPROF test identified two biological zones at Valdelagrana and Cortadura and Lances, and three biological zones at Hoyo, Terrón, Bota, Mazagón, Costa Ballena, Levante, Barrosa and Zahara. Exceptionally, Matalascañas did not show a clear zonation pattern. For this analysis, the sampling levels where no species were present were removed.

A global zonation pattern of the entire set of beaches from the Spanish Gulf of Cadiz coastline was derived from the individual across-shore species distribution; therefore, faunal zones identified at each beach were gathered for a global MDS ordination (Fig. 3). SIMPER analysis performed on this ordination showed a degree of similarity of 40% between all lower zones, where Pontocrates arenarius, Gastrosaccus sanctus and Scolelepis spp. showed the highest percentages of presences (17.8%, 17.2% and 11.0%, respectively). The middle zones showed a similarity of about 30%. Polychaetes of the genus Scolelepis (37.70%), the isopod Eurydice affinis (26.40%), the amphipod Haustorius arenarius (11.56%), and nemerteans (9.95%) highlighted the similarity in faunal composition between all middle zones. Finally, the upper zones showed a 20% similarity and the typifying species were the air-breathing amphipod Talitrus saltator (56.7%) and coleopterans (37%).

Biologically, density values decreased from the lower to the upper zone. In the lower and middle zones the most abundant taxa were crustaceans and polychaetes, while in the upper zones, besides crustaceans, insects were dominant (Fig. 4).

Environmental variables significantly related to fauna variations, tested by Monte Carlo permutation test, were elevation (p=0.002), sand moisture (p=0.001), organic matter content (p=0.015) and grain size (p=0.001). However, these predictor variables were not strongly correlated (r²<0.5). The Monte Carlo test for the set of environmental variables was significant for both axes (p=0.001) and for eigenvalues (p=0.001) showing a significant relationship between biological data and predictor environmental variables.

The CCA results show that the total variation of data was 2.49 (inertia), while the total variation explained was 0.802 (the sum of all canonical eigenvalues). Pearson species-environmental correlations were relatively high: 0.93 for Axis 1 and 0.82 for Axis 2. The first axis (accounting for 66% of the total explained variation) correlated positively with elevation (0.745) and negatively with sand moisture (–0.887) and organic matter content (–0.465). The second axis (accounting for 0% of the explained variation) correlated mainly with medium-sized grain (0.806).

The ordination diagram of CCA (Fig. 5) showed a gradient of zones (lower, middle, and upper) marked mainly by the first axis and showed that crustaceans (Bathyporeia pelagica, Eurydice affinis, E. pulchra, Gastrosaccus sanctus, G. spinifer, Haustorius arenarius and Pontocrates arenarius) and Polychaeta (Scolelepis spp.) responded positively to sand moisture and organic matter content but negatively to elevation, increasing their density to the left along the first axis. Coleopterans and Talitrus saltator exhibited the opposite pattern. Density of nemerteans was the least explained by these environmental variables.

The nemerteans P. arenarius and G. sanctus also responded positively to medium-coarse grain size, while the density of Bathyporeia pelagica, Donax trunculus and Coleoptera sp1 was more influenced by fine grain size, due to their distribution along the second axis.

Table 2 provides a summary of 12 cases of beach usage by means of CI and RI scores and the ensuing management strategy is illustrated in Figure 6. Most beaches are in the multiple-use category, indicating that they have high conservation and recreation values. Only Matalascañas and Valdelagrana are in the intensive recreation category, with high RI values (7 and 9, respectively).

Sandy beaches are extremely dynamic ecosystems with hostile conditions for life, so the number of taxa adapted to living under these conditions is expected to be low compared with other coastal systems. However, the macroinfauna biodiversity found in sandy beaches of the Spanish coast of the Gulf of Cadiz were relatively high (from 4 to 33 species). This value is similar to those reported in nearby latitudes, such as the 9 to 31 species found in northern Spain (Rodil et al. 2006), the 27 to 30 species found on the western coast of Portugal (Gonçalves et al. 2009), and the 4 to 24 species found on Moroccan sandy beaches (Bayed 2003).

The beaches showed a wide range of morphodynamic types and a general increasing trend in species richness from reflective to dissipative beaches was observed, in agreement with McLachlan et al. (1993). Bota beach showed the highest species richness and abundance. This beach, located near the mouth of the Piedras River, where the influence of wave action is lower, is one of the most sheltered of the entire set of beaches studied, as is also reflected in the relative tidal range, which showed high values on this sandy beach. The high richness and abundance values found on Bota beach supports the general trend of biotic variables to increase with exposure, as reported by other authors (Dexter 1983, Jaramillo and MacLachlan 1993, Rodil et al. 2007). Abundance and richness of macroinfauna are directly related to food supply (Rodil et al. 2012), so it is also possible that Piedras river mouth increases food availability, enabling the establishment and development of more species.

Crustaceans, polychaetes and molluscs were usually dominant among the macroinfauna of sandy beaches (McLachlan and Brown 2006). In our study, amphipod and isopod crustaceans and spionid polychaetes were the most abundant and included diverse taxa. In fact, 74% of all individuals collected belong to six species of these groups: Bathyporeia pelagica, Haustorius arenarius, Pontocrates arenarius, Siphonoecetes sabatieri, Eurydice affinis and Scolelepis spp. The dominant and most frequent species occurring on every beach studied were polychaetes of the genus Scolelepis. This genus has a wide geographical distribution (Souza and Borzone 2000) and is also very abundant on many Atlantic beaches (Papageorgiou et al. 2006, Degraer et al. 2003, Barros et al. 2001).

Little importance is given to nemerteans, which are normally not considered typical taxa on sandy beaches, due to residual contributions that they exhibit, although this taxon is considered a useful bio-indicator (Herrera-Bachiller et al. 2008). For example, the nemerteans Lienus acutifrons and Psammamphiporus elongatus were affected by the “Prestige” oil spill on Galician sandy beaches, so they could be a good indicators of beach quality (Herrera-Bachiller et al. 2008, Puerta et al. 2010). On sandy beaches of southwestern Spain, Nemertean abundance was similar to that of molluscs, showing high occurrences (67% of the sampled beaches), highlighting the importance of nemerteans at these latitudes. Similarly, Talitrus saltator was frequently found on the beaches studied. This ‘sand-hopper’ is recognized as a good biomonitor of trace-metal pollution and the effect of human trampling (Ugolini et al. 2008).

The zonation of macroinfauna, defined as the distribution of species throughout the intertidal zone on sandy beaches, has been undertaken around the world (Defeo et al. 1992, Jaramillo et al. 2000, Schlacher and Thompson 2013). Macroinfauna across-shore distribution is highly variable, ranging from one to five biological zones, although three are most common (see Schlacher and Thompson 2013). In the current study, 67% of total beaches showed three distinct zones and 25% showed two. This pattern was different to that found on beaches in northern Spain, where macroinfauna was divided into two main zones (Rodil et al. 2006), and on Moroccan sandy beaches, where four zones were identified according to Salvat’s scheme (Salvat 1964, Bayed 2003).

Jaramillo et al. (1993) determined that intermediate and dissipative beaches include three faunal zones, whereas the reflective beaches had only two. Along the Spanish coast of the Gulf of Cadiz, this pattern was not found. In fact, the more dissipative beaches showed two biological zones, while beaches closest to the reflective state (Hoyo and Mazagón) had three zones. In general, the number of zones alternated, independently of Dean’s parameter. Thus, no clear evidence was found to support the contention that the number of zones is closely related to morphodynamics. These results corroborate the conclusion drawn by Schlacher and Thompson (2013), who detected no significant correlation between habitat metric (habitat dimensions, sediment properties, and morphodynamic state) and the number of fauna zones.

Although the number of biological zones varied among beaches, a common zonation pattern was established for the entire set of beaches studied. This was done in order to characterize the most typical species inhabiting each zone. The general pattern showed three biological zones according to the zonation scheme of Dahl (1952). In general, the supralittoral zone was typified by air-breathing amphipods (Talitrus saltator) and coleopterans. The middle zone was dominated by true intertidal species, such as Haustoriidae amphipods (Haustorius arenarius), Cirolanidae isopods (Eurydice affinis), Spionidae polychaetes (Scolelepis spp.) and nemerteans. The lower or sublittoral zone was typified by amphipods belonging to the Pontoporeiidae family: mysids, and spionid polychaetes. The distribution of the species in each zone is in agreement with findings on other nearby temperate sandy beaches, such as those along the northern coasts of Spain, Tunisia, and Morocco (Rodil et al. 2006, Pérez-Domingo et al. 2008, Bayed 2003).

Diversity and densities of individuals increase towards the lower zones. This is a general feature found in numerous studies of sandy beaches worldwide (Jaramillo et al. 1993, Gonçalves et al. 2009, Rodil et al. 2006). Some authors have determined that this pattern could be due to a reflection of the high sub-tidal diversity and short periods of air exposure, thus allowing more species to inhabit zones closest to the seawater (Degraer et al. 1999, Aerts et al. 2004). The high abundance found in the lower areas of all the beaches studied evidences how important these environments are as potential sources of food for other predatory species (fish and birds).

Distribution of macroinfauna is related to the tolerance of these communities to environmental variables (McLachlan and Brown 2006). Although the relationship between species and the environment could change with the scale of study (Rodil et al. 2012), abiotic predictor variables at the local scale were examined.

Beach slope and grain size have been identified as main factors controlling macroinfauna distribution throughout the intertidal zone (McLachlan et al. 1993, Jaramillo and MacLachlan 1993). Results from CCA analysis showed that sand moisture and organic matter content, in addition to elevation and grain size, were the main environmental variables controlling the macroinfauna distribution across the shore on sandy beaches of the Gulf of Cadiz coast.

Lower and middle zones showed an internal gradient influenced mainly by average grain size. Species inhabiting these zones were Pontocrates arenarius, Siphonoecetes sabatieri and nemerteans, closely related to coarse grain size, while Donax trunculus and Bathyporeia pelagica were related to fine grain size.

The most abundant species in upper zones, the talitrid amphipod Talitrus saltator and coleopterans, were correlated positively with elevation but negatively with sand moisture and organic matter content. Grain size was not a good explanatory variable for these species. In fact, Ugolini et al. (2008) found no relationship between sand-hopper abundance and sand-grain size. Although these species showed significant relationships with abiotic variables, other factors not taken into account could affect their distribution. For example, it has been reported that stranded material (e.g. macrophytes and macroalgae) provides a physical structure that can be used as a shelter or breeding site and as a food source by supralittoral arthropods (Olabarria et al. 2007) and the age of these deposits plays a significant role in the structure of upper-shore assemblages (Ruiz-Delgado et al. 2014).

The conservation and recreation potential indices developed by McLachlan et al. (2013) are useful tools for determining the most suitable beach management approach, providing a set of management principles that can easily be followed, even with limited resources. In these cases, mapping plays a key role for the collection of information and has been identified as a principal step towards management and conservation policies (Rolet et al. 2015). In addition, the implementation of these indices provides the broadest characterization of sandy beaches, not only based on physical features, since sandy beaches are also social and ecological environments.

In this work, 12 case studies illustrate the application of the CI and RI indices. Although these indices include five categories, our beaches are grouped mainly into two: intensive recreation and multiple use. Most beaches have shown scores above 5 on the CI, so they have a reasonable potential for conservation. For example, Terrón, Mazagón and Levante beaches, categorized by these indices as of multiple use, have a high conservation potential. Coastal management policies should therefore be designed to promote measures for the conservation of these environments. On the other hand, only Matalascañas and Valdelagrana, have a high recreation potential, which can be exploited with appropriate management strategies or the possible re-evaluation of existing measures that reduce the CI value, for example.

In conclusion, the beaches on the Spanish coastline along the Gulf of Cadiz are characterized by high biodiversity, including major bio-indicator species, and by a clear zonation of macroinfauna. The overall distribution pattern involves three biological zones: the supralittoral zone, characterized by air-breathing amphipods and coleopterans; the middle zone, characterized by Haustoriidae amphipods, Cirolanidae isopods, Spionidae polychaetes and nemerteans; and the sublittoral zone, characterized by amphipods belonging to the Pontoporeiidae family, mysids and spionid polychaetes. The across-shore distribution of the macroinfauna is influenced primarily by sand moisture, organic matter content, elevation and grain size. Other factors such as wrack deposit and organic inputs from rivers and estuaries could influence the abundance and distribution of macroinfauna inhabiting sandy beaches. Therefore, future studies are needed to elucidate whether the presence of stranded material could affect the global zonation patterns on sandy beaches. The recreation and conservation indices showed that the studied beaches have multiple functionalities, given their social, physical, morphodynamic and biotic features, although many of them have a high potential for conservation that needs to be managed.