Abstract

This study assesses the impacts of wave action and freshwater outflow on soft-bottom benthic macrofauna spatial distribution and temporal stability along the highly exposed French Basque coast. Sediment characteristics and macrofauna abundance have been seasonally investigated during two years for nine stations located at the same (6 m) depth and spread over three subtidal sites showing distinct exposure levels. Wave climate has been determined through an operational numerical model. A total of 121 taxa were recorded, gathered in three main faunal assemblages, as revealed by classification and ordination methods. Non-parametric multivariate multiple regression (distance-based linear model) showed that the variations in macrofaunal distribution can be explained by hydrodynamic conditions. Wave exposure strongly linked to estuarine inputs were the most relevant abiotic factors influencing distributional patterns and functional structure as described by biological traits analysis. Despite the influence of these abiotic variables affecting sedimentary dynamics, seasonal stability was observed in macrobenthic assemblages composition suggesting an ability to recover from natural disturbances such as (e.g.) winter storms. In this way, these results provide baseline knowledge for future ecosystem and resource management in shallow subtidal areas strongly exposed to swell and freshwater outflow where soft-bottom macrozoobenthic communities are less frequently studied.

1. Introduction

With the European Directives (Water Framework and Marine Strategy) adoption, ambitious objectives for the conservation and the restoration of the state of water bodies have been fixed. In order to establish the ecological quality of European coastal and estuarine waters, the Water Framework Directive emphasizes the importance of biological indicators, such as phytoplankton, macroalgae, benthos and fishes [1, 2]. More particularly, macrobenthic fauna is an important component of marine ecological systems being involved in nutrients cycling, pollutant metabolism and constituting a food source for higher trophic levels [3, 4]. Macrobenthos are also known to be an effective indicator. Due to their longevity, sedentary nature, different tolerances to environmental stress, intermediate trophic level positions, and close association with the substrate, benthic macrofauna integrates effects of environmental variations and provides a relatively clear signal, susceptible to detect a disturbance on the ecosystem [59].

Assessing variability in biodiversity along environmental gradients and identifying factors responsible for spatial patterns in macrofaunal assemblages are a central theme in marine benthic ecology [10, 11]. Since Petersen’s study [12], studies conducted worldwide consider sediment characteristics as a major explanatory factor of benthic macrofauna distribution [1318]. Moreover, at broader spatial scales, other natural environmental factors such as the hydrodynamic conditions (including freshwater outflow and wave action) seems to directly or indirectly control the presence and abundance of macrobenthos [19, 20]. Water movements through transport of sediment and organic material, strongly affect soft-bottom community by smothering immobile forms or forcing mobile forms to migrate, altering grain size distribution, impacting light penetration and primary productivity, alternating episodic erosion and deposition processes [2123]. These natural physical disturbances impact the structure of macrobenthic communities and also their functional responses to important ecosystem processes, such as resource usage (nutrients recycling), feeding interactions (trophic structure), habitat building organisms (ecosystem engineering, biogenic structures), bioaccumulation (body size, growth rate, longevity) and sediment properties (tube-building, burrowing activities) [2428]. Therefore, hydrodynamic conditions, as natural physical disturbances, determine the colonization of a given habitat [22], they alter interactions between individuals [29] and consequently shape macrobenthic community [21, 23, 28, 30, 31].

Acting within a different scale of time, from daily to seasonal variations, natural physical factors such as hydrodynamic conditions and their associated sedimentological processes may cause both press and pulse types of disturbance [32]. Press-disturbance process causes troubles by acting over a prolonged period of time that is intolerable to benthos.The intermediate frequency and persistence of the disturbance pressure are higher than the endurance and rate of ecological succession of the biota. On the other hand, pulse process causes a disturbance by exceeding a threshold above which benthos are unable to remain attached to the seabed or are buried under rapidly deposited sediment. On the continental shelf, press-type disturbances include the sustained period of elevated turbidity that follows a storm or a flood event [3335], and pulse-type processes include the mobilization of bottom sediments by storms [36, 37]. Ecosystems with strong physical constraints (i.e. estuaries, shallow habitats) are characterized by low diversity (e.g. species richness) and species with an inherent ecological plasticity allowing them to sustain their domain of stability when facing external disturbances [3841]. Thus, they are also expected to recover quickly following a disturbance than communities in more stable environments such as deep habitats [42].

The Basque country is exposed to strong swell because of its location in the innermost part of the Bay of Biscay and the small width of the continental shelf. This coast is additionally bordered by the estuary of Adour in the North and many small mountainous coastal rivers in the South contributing to a huge proportion of the sediment fluxes into the Bay of Biscay [43]. Numerous studies of soft-bottom macrobenthic community have been conducted along the estuaries and the continental shelf of the Spanish coast [4453]. These communities have been surveyed less frequently for coastal sandy beaches [54]. In these nearshore areas, rocky substrates are indeed more documented than soft-substratum [55]. Regarding the high level of anthropogenic threats and disturbance within these ecotonal zones, and the requirements of EU directives, studies such as the current one are essential to provide baseline knowledge that can be enhanced for sustainable management of these areas [56]. Indeed, the major constrains to implement conservation strategies in marine ecosystems are the general lack of baseline data prior to impacts and the substantial gaps in the current knowledge of natural patterns of variability of their assemblages, which are intrinsically variable [20, 57, 58]. Moreover, assessing the main natural environmental factors that shape spatial patterns of macrobenthic assemblages and their temporal stability will help to discriminate between natural and anthropogenic changes [11, 59] and to appreciate their resilience in a context of applied ecological research expansion on the macrobenthic communities residing in coastal and estuarine areas [58, 6062]. Therefore, the aims of this study were to characterize the ecotonal macrobenthic nearshore soft-bottom communities in an area strongly exposed to wave and estuarine inputs, to determine natural environmental factors that shape spatial patterns and discuss their temporal stability. To achieve these aims, first spatial patterns of sediment features, wave climate, estuarine inputs and macrobenthic assemblages were described within three sites located at the same depth and showing distinct hydrodynamic conditions (exposed coast: Anglet coast, semi-enclosed bay: Saint-Jean-de-Luz Bay, opened bay: Hendaye Bay). Then the relationship between spatial distribution patterns of macrobenthos and those of environmental factors were investigated using multivariate statistical approaches. This will allow identification of useful predictor variables and relating them to biological traits along shallow strongly exposed sandy subtidal coasts, such as the French Basque coast. Finally, the temporal patterns of macrobenthic assemblages were evaluated through time for the whole nearshore area from 2014 to 2016.

2. Material and Methods

2.1. Study Area

This study was conducted along the French Basque coast, located at the south western part of the French atlantic coast (Figure 1). In this area, tide is semi-diurnal and mesotidal, ranging from 1.85 to 3.85 m. Waves predominantly come from the West-North West direction with a 10-s peak period and a 2-m average significant wave height [63].

The studied area is subjected to the influence of three main rivers which are, from the north to the south: the Adour River, the Nivelle River and the Bidassoa River. The Nivelle and Bidassoa Rivers, subjected to flash floods, are regarded as torrential rivers considering their watershed slopes [43]. By contrast the Adour river is characterized by a relatively flatter and larger watershed. In terms of sediment discharge, the Nivelle and Bidassoa solid flow are ten times less important than Adour sediment discharge [64, 65].

2.2. Sampling Procedure and Laboratory Analyses

To assess how wave conditions impact subtidal soft-bottom macrozoobenthic communities, three sites (Anglet coast, Saint Jean de Luz Bay, Hendaye Bay, see Figure 1) subjected to different estuarine inputs and wave exposure levels were seasonally investigated between August 2014 and June 2016. At each of the four seasons (August, December, March and June) the three locations were sampled for a total of eight sampling campaigns per station during the two years studied (sampling dates available in Supplementary materials Table S1). The nine sampling stations (three per sites) were located at the same depth (−6 m chart datum).

To assess the macrobenthic assemblages, three sediment samples were collected using a van Veen grab (0.1 m2). Grab contents were sieved through a 1 mm mesh size. Material retained on the sieve was directly fixed in ethanol (99.9%) for later identification to the lowest taxonomic level (predominantly species) and enumeration in the laboratory. The World Register of Marine Species [66] was used to check and harmonize species names.

For the sediment analysis, a very small sub-sample of each collected grab was used for the determination of both sediments organic matter content and sediments grain size analyses.

2.3. Environmental Data
2.3.1. Wave Climate

Wave climate was determined for each station and between each field campaign from a SWAN operational model developed within the European project Littoral, Ocean and Rivers of Euskadi-Aquitaine (LOREA). Detailed model setup and validation results were further described in Dugor et al. [67]. The model boundaries were forced by HOMERE sea-states hindcast database, based on WAVEWATCH III model. Wind data were provided by the ECMWF (European Center for Medium-Range Weather Forecasts). A nesting strategy allowed making the transition between offshore and coastal models over three successional grids: a regional grid, an intermediate grid and finally three local grids with a 20 m resolution around studied sites. Four wave parameters were obtained in order to describe wave climate: mean significant wave height (Hsmean), maximum significant wave height (Hsmax), mean bottom orbital velocity (Ubrmean) and maximum bottom orbital velocity (Ubrmax). The wave climate characterization was carried out for the period preceding each sampling campaign.

The three-dimensional ECOMARS model [68, 69], as performed by Dutertre et al. [11] at a larger scale, was not used in this study. Resolution grid (3 km) was not adapted to represent hydrological variations within the three sampled sites, where stations were spaced from 400 m in the Saint Jean de Luz Bay to 1.3 km along the Anglet coast.

2.3.2. Sediment Characterization

Data used for sediment characterisation were treated as percentages for each grain size categories determined using a sieve shaker. The following sedimentary fractions were considered based on the classification of Wentworth [70] modified by Folk [71], Folk and Ward [72] and Folk [73]: “GR” gravel and pebble (>2 mm), “VCS” very coarse sand (1–2 mm), “CS” coarse sand (0.5–1 mm), “MS” medium sand (0.25–0.5 mm), “FS” fine sand (0.125–0.25 mm), “VFS” very fine sand (0.063–0.125 mm) and “F” silt & clay (<0.063 mm). The diameter corresponding to the average grain size of sediment particles (D50) and the sorting index (So, [74]) were calculated using a MATLAB routine for each station and each field campaign. D50 was expressed in the phi (φ) scale originally developed by Krumbein [79] in order to simplify statistical analyses. Organic matter (“OM”) content was estimated by loss of ignition (450°C, 6H) and was also treated as percentage of sediments weight.

2.3.3. Estuarine Inputs

The study area is influenced by three main rivers which are, from the north to the south: the Adour River, the Nivelle River and the Bidassoa River. To take into account estuarine influence, mean and maximal river discharge (respectively “Qmean” and “Qmax”) were retrieved between each field campaign from the French water information system database (http://www.hydro.eaufrance.fr/) and from the Confederación Hidrográfica del Cantábrico (http://www.chcantabrico.es/). Distance between each sampling station and the river mouth (“DistMouth”) was also determined as a proxy of salinity level and freshwater influence (with reference to fresh water input; Table 1).

2.4. Data Analysis
2.4.1. Environmental Variables

Analyses of variance (ANOVA) were performed using R Software to test for difference in environmental variables among localities. These analyses were based on a one-way model, including locality as fixed factor with three levels. Assumptions of data normality and homogeneity of variances were previously assessed using Shapiro-Wilk and Levene’ test, respectively. Whenever ANOVA’s assumptions were not met, a non-parametric one-way analysis of variance was performed (Kruskal Wallis’ H test). Post-hoc pairwise multiple comparisons were performed using the Tukey test whenever ANOVA showed significant differences () and the Nemenyi test was used following nonparametric one-way analysis of variance.

2.4.2. Macrobenthic Communities

(1) Spatial Distribution of Soft-bottom Communities.

The structure of the macrobenthic community was investigated using multivariate techniques provided by PRIMER software [75]. Original data consisted of “stations × species” matrix which was obtained after removing rare species. Species were considered as rare when they only appear in a single station and with a contribution to the station total abundance lower than 5%. Abundance data were -transformed prior to analysis.

Similarity relationships between stations of all biotic data results (eight sampling dates) were determined using the Bray–Curtis coefficient [76]. The objective was to assess the spatial distributions of soft-bottom communities taking into account temporal variability of the macrobenthic community at each sampling station. Therefore a similarity matrix was first computed for all the stations-dates. This matrix was then used to compute the matrix of distance among sampling stations centroids using the “distance among-centroids” routine in PRIMER. This procedure allows for comparing the different sampling stations (the centroids) while integrating each stations variability obtained through the different seasonal sampling campaigns (each data point that were used to determine the centroids). The matrix of distance among centroids obtained was used to perform a hierarchical cluster analysis using group-average clustering (in accordance with Legendre and Legendre [77]) in order to identify groups of stations displaying similar fauna. A Principal Coordinate Analysis (PCO) was also performed on the stations centroids to show the defined groups in a two dimensional space [78]. Each benthic assemblage, resulting from the multivariate analyses, was then characterized by its species richness (S), density of individuals (N), Shannon’s diversity index (H’, log), and Pielou’s evenness index (J’). A similarity percentages (SIMPER) analysis was finally used to determine contributions of each specie to the Bray-Curtis similarity within each of the groups.

Relations between environmental variables and benthic community distribution was assessed using distance-based linear models (DISTLM) which consists of partitioning variability in the dissimilarity matrix according to environmental variables as predictors [78]. Ten environmental variables were taken into account: wave climate through mean significant wave height, maximum significant wave height, mean bottom orbital velocity and maximum bottom orbital velocity; estuarine inputs through mean and maximal river discharges and distance to the river mouth and, finally, sediments variables as D50, sorting index and organic content. These relations were assessed using the PERMANOVA + add-on [78] of the PRIMER software. Prior to these analyses, a selection of variables was performed by selecting among the variables displaying high level of Spearman rank correlation coefficient (≥0.9, disregarding the sign of the coefficient).

RLQ analyses were performed to relate environmental variables with a significant influence on macrofauna distribution (coming from the DISTLM results) to biological traits [80, 81]. This method requires the generation of three different data tables: the table gathering information of the significant environmental variables per site and sampling dates according to the DISTLM results, the table with the abundance of each specie in every samples for different dates and the table composed of biological trait data [27, 80, 82]. Information about functional traits was compiled by gathering information from several literature sources: species identification guides, research papers and web database. The main limitation of Biological traits analysis (BTA) is the occurrence of gaps in the knowledge of some species’ biology [27, 83]. To minimize this aspect, the list of species in the BTA analysis was reduced without causing the loss of integrity in analysis. Thus, from a total number of 121 taxa, 101 species were considered which contributed to 95.5% of the total abundance observed. All species characteristics of faunal assemblages determined by SIMPER procedure showed biological traits description. Each trait was divided into a maximum of five modalities representing different categories of a trait displayed by the considered organisms (Table 2). The open source software R and the Ade4 package [84] were used to perform the RLQ analysis.

(2) Temporal Stability of the Observed Faunal Patterns.

In order to measure the temporal variation in faunal assemblages at the scale of the whole nearshore area, the abundances of individual taxa were averaged across the sampling stations for each faunal assemblage at each sampling time (8 field campaigns). For each assemblage, patterns of dissimilarity through time were visualized using a Principal Coordinate Analysis (PCO) of the assemblage × time centroids. The temporal variability (dispersion) of each assemblage was quantified as the average Bray-Curtis dissimilarity among time points. These dispersions were formally compared among the 3 assemblages using a permutation test of dispersion with 9999 permutations (PERMDISP, see [85]). This approach directly compared temporal variation in the community structure of whole assemblages of the nearshore area. In addition, to compare the station-level temporal variation among assemblages, the average and the standard error of temporal variation calculated from the stations were plotted for each faunal assemblage.

3. Results

3.1. Environmental Variables
3.1.1. Wave Climate

Data on wave regime are compiled in Figure 2. The one-way ANOVA and multiple comparisons for the wave exposure level are available in Supplementary material (Tables S2a, b). Significant differences among the 3 localities were observed for the four wave parameters. Wave exposure was higher along the exposed Anglet coast (range 0.78–1.96 m), intermediate in the opened Hendaye Bay (range 0.47–1.33 m) and relatively low in the semi-enclosed Saint Jean de Luz Bay (range 0.28–0.78 m).

There was no significant difference of wave climate among the sampled stations within each of the three sites (Tables S3, S4 and S5).

3.1.2. Sediment Features

Data on sediment composition and organic matter content (OM) were compiled for each stations during the eight field campaigns (Table 3). The one-way ANOVA and multiple comparisons tests for sediment parameters among the three localities are available in Supplementary material (Tables S6a, b). Significant differences were observed for all sedimentary parameters excepted for gravel content. Post hoc analyses indicated that grain size along the Anglet coast was significantly coarser and globally better sorted than in the remaining localities. Additionally the Anglet coasts’ stations showed a significantly lower concentration of OM in sediments.

Along the Anglet coast, the northern station (A12) consisted of clean, medium to fine sands with gravels and coarse sands; whereas the two other stations (A13, A14) consisted of clean, fine sands. The mean level of OM was lower than 1% at all stations in the Anglet coast (Table 3). In the Saint-Jean-de-Luz Bay, the eastern station sediments (N5) consisted in slightly muddy, heterogeneous sand, clearly coarser than the 2 other stations which consisted of muddy, fine sands. Organic matter content was relatively high, with an average up to 5% (Table 3). In the Hendaye Bay, the station located closer to the river mouth (B11) showed the highest proportion of silt and clay fraction (up to 10%) close to the one observed in the western part of the Saint Jean de Luz Bay. This station contained also the highest level of OM (Table 3) compared to the two others sampled in this site. These differences among sedimentary features were significant within the Saint Jean de Luz and Hendaye sampled stations (Tables S7S9 in Supplementary material).

3.1.3. Estuarine Inputs

Heterogeneity appeared among the three localities regarding freshwater influence (Supplementary material: Table S10). The Adour river displayed a much more important mean daily river discharge of 311 m3·s−1 compared to less than 120 m3·s−1 for the Bidassoa and less than 30 m3·s−1 for the Nivelle (Figure 3). Stations located along the Anglet coast are however farther from the river mouth (more than 1.8 km) and located along an open coast compared to the stations of the St-Jean-de-Luz and Hendaye Bays which are probably more directly impacted by fresh water inputs because they are closer to the river mouth and located in embayments.

3.2. Macrobenthic Community Distribution

A total of 121 taxa were recorded at the nine sampling stations during the eight seasonal campaigns. Crustaceans and polychaetes were the most diverse groups with respectively, 45 (37%) and 40 species (33%). Molluscs included 26 species (22%) and echinoderms eight species (7%). Other species (1%) belonged to Nemertea. The most frequently observed species were the crustacean Diogenes pugilator and the polychaete Nephtys cirrosa respectively found in 75% and 55% of the records.

The dendrogram produced by the hierarchical agglomerative clustering distinguished three main groups of stations corresponding to three main different species assemblages (Figure 4): the first dichotomy of the dendrogram, separated species assemblage A from the other two assemblages (B and C), which displayed a higher level of similarity (Figure 4).

Species assemblage A included only samples from the exposed Anglet coast (northern part of the investigated area, Figures 1 and 4). The top two contributive species (Polychaetes Scolelepis spp., Nephtys cirrosa) are typical species of exposed sandy bottoms. This assemblage displayed the lowest average species number and density compared to the two other assemblages identified (Table 4).

Species assemblage B included stations located in the western part of the Saint Jean de Luz Bay and in the middle part of the Hendaye Bay (Figures 1 and 4). This muddy sand assemblage showed the highest species number and fauna density with 25 ± 9 species and 162 ± 211 individuals 0.3 m−2 (Table 4). This assemblage was mainly characterized by Diogenes pugilator together with the molluscs Fabulina fabula, Abra alba, Tritia reticulata and Antalis novemcostata, the Polychaetes Nephtys hombergii and Sigalion mathildae, as well as the sea urchin Echinocaridum cordatum.

Because it shared some contributive species with both previously described assemblages (Diogenes pugilator, E. cordatum, T. reticulata: from assemblage B; N. cirrosa: from assemblage A), assemblage C composition can be considered as intermediate between these previously described assemblages. The level of fauna density and number of species were also intermediate between those measured in the previously described assemblages. Regarding sediment type, this assemblage was retrieved from slightly muddy sands (intermediate between sands characterizing assemblage A and muddy sands characterizing assemblage B, see Table 4) and located at the entrance of bays: the eastern and the western parts of the Hendaye Bay (stations B12 and B10) as well as station N5 located at the eastern part of the Saint-Jean-de-Luz Bay (Figures 1 and 4).

3.3. Environmental Drivers of Macrofaunal Distribution

Maximal significant wave height was, by far, the main variable explaining the highest amount of variance in macrofauna abundance (Table 5). Alone, this parameter explained 47% of the variance, followed by organic matter (12%), D50 (8%), and sorting index (10%). The last three variables taken individually were not statistically significant (). The best model to explain the macrofaunal distribution would include only the wave exposure parameter. It should be noticed that in this dataset, maximal significant wave height (Hsmax) was positively correlated to the three other hydrodynamic variables (Hsmean: 1.00, Ubrmean: 0.98, Ubrmax: 0.98) and the estuarine inputs parameters (DistMouth: 0.87, Qmean: 0.97 and Qmax: 0.96).

3.4. Relation between Traits Modalities and Environmental Variables with a Significant Influence on the Faunal Distribution

To relate significant environmental variables explaining the faunal distribution to biological traits, a RLQ analysis was performed.

This analysis identified the high hydrodynamic area. Strong associations were observed between the positive part of RLQ axis 1 and the maximal significant wave height (Hsmax) and between the negative part of RLQ axis 2 and the maximal river discharge (Qmax). The corresponding biological traits associated with higher exposure were walking motility, demersal and free-living epifaunal species indifferent to sediments organic matter enrichment (AMBI group II) (Figures 57).

In contrast, the negative part of the RLQ axis 1 and the positive part of the RLQ axis 2, associated with lower exposure level, depicted sessile motility, surface and sub-surface deposit-feeders and tolerance to opportunistic species.

3.5. Temporal Stability of These Faunal Assemblages

Individual benthic assemblages were identifiable as clusters of points having similar symbol and color on the PCO plot (Figure 8). The temporal variation of any individual assemblage is measured and can be seen in two ways in this figure: (i) from the relative spread (dispersion) of time points for each assemblage in the PCO plot and (ii) from the bar graph showing the average Bray-Curtis (BC) dissimilarity among time points for each assemblage.

The individual assemblages globally formed distinct clusters on the PCO plot (Figure 8). The first axis of the PCO corresponds to a clear gradient of exposure among assemblages: the most wave exposed assemblage (assemblage A) was located on the positive part of axis 1 while samples corresponding to the more sheltered assemblage B were located on the negative part of the axis.

Despite an apparent greater degree of dissimilarity among samples collected on different dates for the assemblages located on wave exposed shallow-water and on opened bay (assemblages A and C, respectively) than for the most sheltered assemblage B, no significant difference in dispersion (temporal variability) was observed among faunal assemblages (PERMDISP, , ). Almost all of the assemblages indeed showed average BC dissimilarities through time around 35% (Figure 8, upper bar graph).

4. Discussion

Responses of benthic organisms to environmental stressors are the integrated result of both direct and indirect processes which can be manifested as changes in abundance, diversity and fitness of individuals and communities [86]. Identifying and integrating the effects of natural pressures is an essential challenge for understanding and managing coastal biotic resources [87, 88] particularly when they are subjected to anthropogenic threats. In shallow subtidal, previous studies have focused on analysing patterns of macrobenthic assemblages along salinity or depth gradients [28, 60, 62, 89]. This study is novel in that it characterizes the benthic macrofauna distribution patterns across three nearshore soft-bottom sites, located at the same depth (6 m) and exposed to different hydrodynamic conditions, along the French Basque coast. The results indicate that environmental variables (wave climate, sediment parameters and estuarine inputs) vary significantly among localities which is common sense, with a clear distinction arising from hydrodynamic conditions, which is a new finding. Anglet coast is indeed clearly distinct from the southern sites with the highest hydrodynamic conditions (e.g. wave exposure and river discharge) and the lowest silt and clay and organic matter contents. Distinctions appeared also within localities within the two bays where sedimentary gradients from East to West were observed.

Relationships between macrobenthic abundance and environmental factors are not easy to explain because they differ among areas [58, 90]. No single mechanism explains patterns observed across many different environments [29]. In this study, the differences in spatial distribution of the 121 taxa found among the whole nearshore area could be 47% explained by wave climate. This correlation between spatial distribution of macrobenthos and natural biotic factors is relatively high. Recent previous studies carried out along the Atlantic coast showed similar degree of variation explained by environmental variables. Along the subtidal coastal fringe of South Brittany, Dutertre et al. [11] found a 51% correlation with a combination of 16 natural abiotic factors whereas Carvalho et al. [28], Veiga et al. [58] and Martins et al. [17] showed significant correlations varying from 35% to 66% on the Portuguese continental shelf. In addition to these published results, this study provides therefore a consistent and thorough understanding of the causes of macrofaunal spatial patterns at the scale of shallow strongly exposed sandy subtidal coasts.

One key finding in this study was that maximal significant wave height explains the largest part of the faunal spatial distribution (47%). Highly linked to freshwater outflow and bottom orbital velocity, these hydrodynamic factors appeared as key descriptors for the local distribution of soft-bottom communities as well as determinant for the sedimentation processes and, consequently sediment types. Indeed, the faunal assemblages gathered the stations composed by similar sediment type. Different community structures were therefore observed within each bay. The less exposed stations of the semi-enclosed bay appeared more similar to the muddy sand stations of the open bay (assemblage B). Conversely, the more exposed station presented a community structure similar to the one observed within the sandy stations of the open bay (assemblage C). Such correlation between macrofauna distribution patterns and the hydrodynamic regime had previously been reported elsewhere, e.g. in Portuguese Continental shelf [17, 28] and in South Brittany [11]. Hydrodynamic conditions, broadly defined as the duration of wave-induced sediment remobilization, were also the most relevant factor, explaining the highest percentage of spatial variation in the macrofauna along the southeastern Portuguese coast [28].

Interactions between benthic organisms and abiotic factors, influencing their environment, result in a wide variety of functional adaptations [91]. Therefore, assessing changes in the functional composition of benthic assemblages in space may shed light on the consequences to the ecosystem services resulting from single or multiple disturbance events [27]. In shallow subtidal areas exposed to strong hydrodynamic conditions, physical erosion and suspension of soft sediment favour infauna and active burrowers, which are the dominant biological traits at the scale of the whole nearshore, as observed in the present study and by Allen and Moore [92] and Dutertre et al. [11]. Suspension feeding (SF) strategy was also the main feeding guilds within the whole environmentally-stressed areas. Suspension feeder communities are indeed generally associated to spaces with strong hydrodynamics acting on the seafloor [9396]. This is related to their dependence on higher oxygen concentrations and the need for small re-suspended particles for feeding purposes [93, 9597]. Within the study area, the RLQ analysis results highlighted nonetheless a pattern of change in trait composition from highly energetic zone to the more sheltered ones. The macrobenthic communities of the highest exposed site (Anglet coast: assemblage A) presented the highest relative densities of free-living fauna such as the swimming crab Portumnus latipes. Within lower hydrodynamics areas, surface deposit feeders (SDF), as Abra alba, were more abundant. This trophic group is generally associated with areas with lower hydrodynamic action on the seafloor, as currents limit their feeding and locomotion abilities [94]. SDF were distributed with higher density in the opened and semi-enclosed bays (Hendaye and Saint Jean de Luz Bays) thanks to lower wave exposure. Such findings are highly generalisable as similar conclusions were made elsewhere [97, 98], species changing their trophic strategy in response to flow and food flux conditions [29, 97, 98]. SDF feed directly on newly deposited organic matter, while sub-surface deposit feeders (SSDF) primarily feed on older organic matter [93, 96, 99]. Within the study area, this guild (SSDF), through for example the species Antalis novemcostata, was associated with the lowest wave exposure of the semi-enclosed bay (Saint Jean de Luz Bay). Similar patterns were depicted by Dolbeth et al. [97] along the Atlantic Portuguese coast. An ecological shift toward “AMBI group II” to “AMBI group IV” was also observed with decreasing sediment grain size and wave exposure level. Second-order opportunistic species was represented by the polychaetes Aphelochaeta sp. and Lagis koreni and the bivalve Corbula gibba which were exclusively sampled within the semi-enclosed bay of Saint Jean de Luz. This result reflected that in naturally organically rich area such as bays close to river mouth, communities generally included opportunistic species and taxa that can also be found in anthropogenically-organic rich areas [5, 100]. Therefore, as a general matter, this study added to other published so far indicate that the spatial pattern of functional composition of benthic species can be used to infer the influence from physical disturbance exerted on ecosystems [27].

Estuarine inputs, as a source of suspended particulate matter, and wave exposure, by inducing sediment remobilization, influenced in opposite directions the grain-size distribution and the organic matter contents of sediments [28]. As confirmed by this study, this natural sediment mobility due to hydrodynamic conditions is a factor of importance on nearshore soft-bottom, controlling the spatial distribution of many species [21, 97, 101] by causing both press- and pulse- types of disturbance as defined by Harris [32]. The effects of such physical disturbances may vary with intensity and duration, and give sometimes dramatic damage to benthic communities, followed by recovery. Nonetheless, nearshore area as a whole demonstrated a rather temporal stability in community structure. Similar results were observed in exposed shallow water worldwide [102104]. These results could be explained by two non-mutually exclusive processes: (i) an increasing number of repeated pulse disturbances (e.g. ocean storms, extreme swell regimes) which can gradually move the system closer to a press response with benthic assemblages being less sensitive [105]; or (ii) the resilience of an ecosystem depends on thresholds of intensity and/or prevalence of the disturbance, but also on the characteristics of the species affected [106]. Highly motile or dispersing species will recover from disturbance faster than the ones with opposite traits [107]. Therefore, it may be proposed, as a general pattern, that macrofaunal communities in shallow exposed subtidal area are composed of species showing high affinities to high hydrodynamic conditions. The most exposed species assemblage located along the Anglet coast, especially illustrates this apparent resilience. Consequently, as demonstrated by Dutertre et al. [11] and in the present study, large-scale ecosystem-based approach (i.e. site to site comparisons) improves the understanding of the relationship between species distribution and environment, and provides a consistent baseline compatible with management concerns and the detection of spatial and temporal changes.

5. Conclusions

Congruent with other published studies, this contribution supports a priori common sense hypothesis that benthic organisms exhibit, as distinct responses to different levels of disturbance [27]. As expected, this study confirms that distinct hydrodynamic conditions do affect the spatial distribution and the functional structure of macrobenthic fauna as revealed from the extensive study of nine subtidal stations (three per site) seasonally sampled during two years and located at the same depth. Environmental constraints represented by wave exposure level and estuarine inputs appear as the more determinant variables for the benthic fauna in narrow environment including these three nearshore soft-bottom sites. Furthermore, the temporal stability observed at the whole nearshore area scale suggests an ability to recover from a natural disturbance.

Therefore, assessing changes in the distribution of benthic assemblages and the functional diversity of these species in space and time may shed light on the consequences from single or multiple disturbance events on the resulting ecosystems. On an applied perspective, knowledges about these biological processes are useful for coastal management topic, allowing to distinguish between natural from anthropogenic variability of macrobenthic compartment in response to disturbances. In shallow, strongly exposed sandy beaches where coastal erosion is a central concern, such baseline information may be particularly determinant to sustainably manage dredged sand dumping and shoreface nourishment.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This work was partly funded by the Region of Nouvelle-Aquitaine and the Chambre de Commerce et d’Industrie de Bayonne Pays Basque.

Acknowledgments

The authors gratefully would like to thank the Department of Pyrenees Atlantiques and the Chambre de Commerce et d’Industrie de Bayonne Pays Basque which lend their boats and made their staff available to carry the field campaigns.

Supplementary Materials

Table S1: Dates of sampling campaigns and their abbreviations. Table S2a: Kruskal–Wallis (One-way analysis of variance) for the mean significant wave height (Hsmean) with post-hoc pairwise multiple comparisons (Nemenyi test) among the three studied localities. Table S2b: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the wave parameters among the three studied localities. Table S3: Kruskal–Wallis (One-way analysis of variance) for the wave parameters among the three studied stations of Anglet site. Table S4: Kruskal–Wallis (One-way analysis of variance) for the wave parameters among the three studied stations of St Jean de Luz Bay. Table S5: Kruskal–Wallis (One-way analysis of variance) for the wave parameters among the three studied stations of Hendaye Bay. Table S6a: Kruskal–Wallis (One-way analysis of variance) for the sedimentary parameters with post-hoc pairwise multiple comparisons (Nemenyi test) among the three studied localities. Table S6b: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the sedimentary parameters among the three studied localities. Table S7a: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the sediment parameters among the three studied stations of Anglet site. Table S7b: Kruskal–Wallis (One-way analysis of variance) for the sediment parameters with post-hoc pairwise multiple comparisons (Nemenyi test) among the three studied stations of Anglet site.Table S8a: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the sediment parameters among studied the three studied stations of St Jean de Luz Bay. Table S8b: Kruskal–Wallis (One-way analysis of variance) for the sediment parameters with post-hoc pairwise multiple comparisons (Nemenyi test) among the three studied stations of St Jean de Luz Bay. Table S9a: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the sediment parameters among the three studied stations of Hendaye Bay. Table S9b: Kruskal–Wallis (One-way analysis of variance) for the sediment parameters with post-hoc pairwise multiple comparisons (Nemenyi test) among the three studied stations of Hendaye Bay. Table S10: One-way ANOVA with post-hoc pairwise multiple comparisons (Tukey test) for the estuarine outputs parameters among the three studied localities. (Supplementary Materials)