Soft-bottom fishes and spatial protection: findings from a temperate marine protected area

Study area

Located on the western coast of Portugal, the Luiz Saldanha Marine Park (LSMP) comprises 38 km of coastline characterized by high rocky cliffs alternating with sheltered bays. With an area of 53 km², the park includes many sublittoral habitat types, including hard and soft substrates. Soft bottoms replace the coastal rocky reefs in most of the area at depths greater than 15–20 m. A substantial mud component is found in the sediment at depths greater than 30 m, while at shallower depths the dominant substrate is sand (medium, coarse and gravely sand) (Henriques et al., 2014). One interesting aspect to point out is the wide variety of soft-substrate habitats and benthic communities found in this area, as described by Henriques et al. (2014). The area has two nearby estuaries (Sado and Tejo) and is next to the Arrábida mountain chain, which shelters it from the dominant north winds. This region is also a biogeographic and oceanographic transition zone between warm and cold temperate waters (Henriques, Gonçalves & Almada, 1999; Lima et al., 2007) which, along with the variety of available habitats (Henriques et al., 2014), makes it an important biodiversity hotspot (Henriques, Gonçalves & Almada, 1999; Gonçalves, Henriques & Almada, 2002).

The MPA was created in 1998 as part of the Arrábida Natural Park, but it was only in August 2005 that the management plan started to be implemented. The management regulation main goals are to enhance the conservation of marine species and promote sustainable fisheries management (Gonçalves, Henriques & Almada, 2002). Given the diversity of human activities, including recreational fishing, scuba diving and commercial fisheries, and the presence of two important fishing harbours: the port of Sesimbra, located within the MPA area, and the nearby port of Setúbal, spatial management was crucial. To mitigate conflicts related to multiple stakeholders and their divergent concerns (Vasconcelos et al., 2013), the LSMP regulations were established gradually over a four-year period and full implementation was achieved in August 2009 (implementation: 23rd Aug. 2005–23rd Aug. 2009; Year 1: LSMP area under buffer area protection level; Year 2: two PPAs were implemented; Year 3: two more PPAs were created, sampling started; Year 4: half of the FPA implemented; 23rd Aug. 2009—the FPA is fully implemented). The regulations define different zones and levels of usage (Fig. 1). Trawling, dredging, purse seining and spearfishing are prohibited in the entire area, and only licensed commercial vessels less than 7 m long are allowed. The fully protected area (FPA) (4 km²) excludes all human activities, with the exception of research and monitoring. The four partially protected areas (PPA), totalling 21 km², allow some recreational activities, including recreational diving, and commercial fishing for cephalopods with octopus traps and jigs, although fishing is not allowed at distances less than 200 m from the shore. The three buffer areas (BA), totalling 28 km², have the lowest level of protection. In addition to octopus traps and jigs (here permitted close to the shore), longlines, trammel nets, gill nets and recreational fishing are allowed (nets at a distance beyond 463 m from the shoreline).

Data collection

Data were collected between November 2007 and May 2012 over a five years period: the first two years during the management plan implementation (period 1—implementation period: August 2007–August 2009) and three years after the MPA was fully implemented (period 2—after MPA full implementation: August 2009–August 2012). Sampling of the three protection levels occurred seasonally, each spring and autumn. Ten sampling campaigns (experimental fishing) and 183 trammel net sets were undertaken (Fig. 1) (total number of samples: 179; four samples were excluded from the analysis due to problems with the fishing procedure). Sampling area covered the FPA (n = 63), two adjacent PPAs (n = 61) and the larger BA (n = 59). To avoid spatial heterogeneity, areas in the eastern and western edges of the MPA were not included (west: proximity to cape; east: proximity to nearby estuary). Four sampling campaigns were undertaken during period 1 (n = 57) and six during period 2 (n = 126) (first three campaigns—four days each, 36 samples; following campaigns—seven days each, total of 147 samples; three samples per day). Two substrate types were sampled: sandy bottom, at 10–20 m (n = 112), and muddy bottom, at 30–45 m (n = 71) (samples: sand-BA: 36, PPA: 38, FPA: 38; mud-BA: 23, PPA: 23, FPA: 25). The BA was considered as the control in the spatial analysis of protection effect due to the allowed use of small-scale fishing gear: gill nets, trammel nets, longlines, jigs, pots and traps. The area surrounding the MPA was not chosen as a control because of greater depths and differences in bottom type. Each sample consisted of the catch of one monofilament trammel net (500 m long; 1.6 m high; inner panel stretched mesh of 100 mm; outer panel stretched mesh of 600 mm; soak time 22–24 h). Nets were set after sunrise and hauled about 24 h later. Each fish was identified on board to the species level and total length (TL) was measured to the nearest millimetre. Specimens were released afterwards, following ‘catch and release’ practice (The Portuguese Institute for Nature Conservation and Forests—ICNF approved the field surveys for this study; ICNF reference no 00540 140307). Biomass was estimated using published length-weight relations (Table S1—list of species with length-weight references, environment, maximum size and fisheries category). Catch in number per unit effort (CPUE) and biomass per unit effort (BPUE) were calculated considering the unit effort as 500 m of trammel net and a soak time of 24 h. This provided average values (abundance and biomass) per sample that enabled comparisons of the three protection levels over time.

Data analysis

CPUE expressed as number of fish per 500 m of net (n⋅500 m⁻¹) was considered as an index of relative abundance and used along with biomass (BPUE—kg⋅500 m⁻¹) to evaluate differences between protection levels and trends over time (before and after full MPA implementation). Prior to data analysis, pelagic species (Auxis rochei, Boops boops, Sarda sarda, Sardina pilchardus, Scomber spp., Trachurus spp.) were excluded from the dataset. Pelagic species are expected to be less meaningful in habitat characterization (Fréon & Misund, 1999) and therefore were not analyzed.

We used generalized linear mixed models (GLMM) (Breslow & Clayton, 1993) to test significance of protection level (BA, PPA, FPA), substrate (sand, mud) and period (period 1, period 2) as predictors of abundance and biomass. The interaction of protection level and substrate was included in the models to assess whether the effect of protection level was different in each substrate type. Sampling campaign (n = 10) was included as random variable to deal with correlation between observations and with the variability related with campaigns (Zuur et al., 2009). Models were fitted by maximum likelihood with Laplace approximation. Poisson and Gamma distributions (log link) were used to fit, respectively, abundance (discrete) and biomass (continuous) data. Variable selection was carried out by likelihood ratio tests between nested models and by looking to the deviance explained by each variable, while model validation was carried out with visual inspection of the residuals (Zuur et al., 2009).

To analyse species spatial and temporal patterns, ratios of abundance, biomass and total length per species (for stingrays—Myliobatis aquila—dimension used was disc width) were also estimated and applied in accordance with a Before-After-Control-Impact (BACI) analysis approach (Underwood, 1992; Underwood, 1994). Comparisons were made between the buffer area (BA) and the two higher protection levels combined (PPA and FPA: PA—protected area), and between period 1 (implementation period) and period 2 (after full implementation) (protection level ratios included data from both periods; period ratios included data from the three protection levels). This specific analysis was done for the species with frequency of occurrence higher than 15% (Table S2). The average ratios (abundance, biomass and total length), standard errors (SE) and 95% confidence intervals (CI) were obtained by bootstrapping (9,999 permutations). The 95% confidence intervals were calculated using the adjusted bootstrap percentile method and used to check for significance (Payton, Greenstone & Schenker, 2003). To avoid problems resulting from small sample size, bootstrapping of total length was only performed for species with more than ten individuals per treatment (protection level or period). Hence, for total length comparisons, the species Citharus linguatula, Raja undulata and Solea solea were excluded from the protection level ratios, and C. linguatula, R. undulata and Rostroraja alba were excluded from the period ratios.

Multivariate procedures were applied to examine differences in fish assemblage according to protection level (PL) and time (period and year). The Hellinger transformation was used for abundance and biomass in order to overcome the problem of zero inflation (Legendre & Gallagher, 2001). Permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001) was performed to test significance of protection level, period and year as factors for community differences. Each constraint factor was assessed individually and then under interaction models: ‘protection level × period’ and ‘protection level × year’. Pairwise analyses were also obtained by comparing the protection levels within each period and comparing periods of each protection level.

The Euclidean distance matrix of transformed abundances was used in a principal components analysis (PCA). This technique provides visualization of multivariate data structure. It was used to explore the similarity between samples according to protection level and period. In order to compare samples, centroids per protection level and period, and dispersion ordiellipses (95% confidence interval) according to the weighted average scores of each protection level, were plotted. A second graphical output was also obtained with the species vectors in order to assess species contributions for community differences. Given the similarity of results obtained for abundance and biomass, the multivariate visualizations based on biomass are not presented.

To explore the main factors structuring fish assemblages, multivariate regression tree (MRT) analysis (De’Ath, 2002) was carried out using the Hellinger transformed abundances. The factors ‘protection level’ and ‘period’, and two environmental variables—depth and temperature were included. Substrate type was not included because it is correlated with depth, and priority was given to the variable that incorporated more information (depth). This analysis creates clusters of sites by splitting the data according to explanatory variable threshold values and site similarity. The clusters and their dependence on the environmental data are represented graphically by a tree. Through this procedure, each cluster defines a species assemblage, and its environment type is described by the associated environmental values. MRT was performed given its ability to deal with discrete variables (in our case, protection level and period), non-linearity, and higher-order interactions among explanatory variables (De’Ath, 2002). The criterion for final model selection was minimization of cross-validated relative error (CVRE).

Indicator values (IndVal) (Dufrêne & Legendre, 1997) were calculated for each species in each node and terminal node (leaf). With this index, species are considered as indicators of a certain cluster of sites according to their relative abundance and relative frequency of occurrence within that group. Specifically, the product of relative abundance and relative frequency of occurrence of the species within a group is calculated and then multiplied by 100. The index takes the value of zero if there is no occurrences of the species within a group, increasing to 100 if the species occurs at all sites within a group and does not occur in any other group. A minimum of 10% contribution to the explained variance was required for a species to be considered as discriminant at a particular node and permutations were used to test for significance. Permutations were also used to assign significant discriminant species for each terminal node. A species was considered an indicator of a certain assemblage according to the cluster for which it had the highest IndVal.

All analyses were done using the R statistical software version 2.15.2 (R Core Team, 2014). The package vegan (Oksanen et al., 2014) was used for the multivariate analysis and the MRT was built with the package mvpart (Therneau & Atkinson, 2012).

The soft-bottom fish community and the MPA effect

In general, our study showed spatial heterogeneity in fish assemblages. Fish abundance and biomass differed according to protection level, with higher values in the PPA and FPA in comparison with the BA. Our results also point out to the importance of substrate, which in association with depth, plays a major role structuring the communities. This is common in coastal assemblages (Claudet et al., 2006; Claudet, García-Charton & Lenfant, 2011; Pais et al., 2014). It is worth noting that not all spatial complexity was incorporated by depth, substrate and protection level, as these variables do not explain the uneven distribution of some species. One outcome of this is that in the sandy substrate, the PPA had higher abundances than the FPA, while in the muddy bottom, the FPA had higher values than the PPA. This pattern is different regarding biomass, with the sandy FPA exhibiting values similar to the PPA. S. senegalensis, which was more abundant in the PPA, and R. clavata (large species with significant contribution in biomass), commonly found in the FPA, are among the species that contributed to these trends (next section moves on to discuss species patterns). Overall, fish assemblages were mainly structured by depth and substrate, followed by protection level. With respect to protection level, some differences in abundance and biomass were found between the PPA and FPA, and the BA was noticeably a distinct group with lower values.

Examining community structure, the two soft-bottom habitats showed different patterns according to protection level. In the sandy habitat, the BA assemblage was found to be distinct from the PPA and FPA, while in the muddy substrate, the PPA and BA constituted a group apart from the FPA. Many studies have demonstrated that partial protection generally results in different communities in relation to no-take areas (Denny & Babcock, 2004; Lester & Halpern, 2008; Guidetti et al., 2014), and this was also what we observed for the muddy stratum. Habitat heterogeneity (Henriques et al., 2014) is probably related with these differences. Another possible factor is the occurrence of illegal fishing, which was reported by Cunha et al. (2011) for this area. The MPA design is another possible cause for the similarity between the muddy PPA and the BA. Because the central PPA corresponds to the area where the marine park is narrowest, about one third of this area does not include depths greater than 30 m, and the corresponding section of muddy bottom is found outside the park’s border. This discontinuity of habitat protection is expected to affect longitudinal migrations of individuals and consequently, the increase of exposure to fishing is likely to prevent effective protection benefits (Gaines et al., 2010).

Conversely, the different pattern that we found in the sandy habitat, with the PPA similar to the no-take community (also regarding abundance and biomass), might be reflecting the effectiveness of the restrictions prevailing in the PPA, meaning that the prohibitions of fishing with longlines, gill nets and trammel nets, i.e., static gear with moderate to low selectivity (Erzini et al., 2003; Erzini et al., 2006; Gonçalves et al., 2007), are effective management measures. This supports the idea that partial protection may offer effective benefits for soft-bottom fish species. This is clearly conditional on the implemented restrictions (Sciberras et al., 2013; Horta e Costa et al., 2016), and fishing gear that are generally regarded as artisanal might also have considerable impact (longlines, gillnets, trammel nets). Other factors possibly contributing to the similarity between partially protected areas and the reserve are related to geographical proximity (favoring potential movement of individuals between these areas—spillover) and site-specific effects (habitat heterogeneity; fishing pressure not equally distributed before the MPA establishment).

To surpass the complexity related to habitat differences and other site-specific effects, additional attention was given to temporal analysis. Comparisons between periods showed that although no significant abundance differences were observed, a biomass increase occurred in the PPA and FPA, with higher values after full implementation of the MPA. Fish biomass increase was the first positive effect observed in many marine reserves (García-Charton et al., 2008; Watson et al., 2009; Di Franco et al., 2012), and it was recently detected in the rocky reef fish community of this particular MPA (Horta e Costa et al., 2013). This trend is possibly related to the decrease in fishing mortality and increase in longevity, allowing fish to grow to larger sizes. Fish assemblages were also analyzed per period. Even though the interaction term of protection and time was not significant, dissimilarities between periods were found both in the PPA and FPA, suggesting that assemblages are possibly being influenced by management measures. Optimal design should have included sampling prior to the management plan implementation (Claudet & Guidetti, 2010; Donovan et al., 2016). First samples were collected with some regulations already in place and protection effect may have been larger than what our data suggests.

The MPA effect at the species level

Overall, a wide variety of species trends and responses to spatial closure was found. Examining the three abundance and biomass response categories (positive, neutral, negative), two subgroups might be considered within the category with positive signs: species with a clear positive response (increase both in abundance and biomass) and species with only one positive indicator (increase only in abundance or biomass). The later includes M. surmuletus, which increased only in abundance and was not significantly more abundant in the higher protection levels, and R. undulata, which showed a biomass increase but no significant abundance trend. The lack of consistency in these results is likely related to the overall low abundance of these species (R. undulata is currently listed as endangered; IUCN, 2017), and a longer period of protection could potentially reveal more information of their trends. In this sense, they are similar to some neutral species, as some showed size increases (M. merluccius, B. podas and M. mustelus), others were more abundant in the PPA and FPA (T. torpedo, M. aquila, C. lastoviza, H. didactylus and S. cabrilla), and others were both more abundant in the higher protection levels and showed increases in size (S. solea).

The group with consistent positive indicators (abundance, biomass and size) includes two rays; R. alba and R. clavata; one gurnard: C. lucerna; and two flatfishes: M. azevia and C. linguatula. For R. alba, this study is the first to confirm spatial protection benefits. This is of particular importance given that this species is classified as endangered (IUCN, 2017). Spatial management seems to be a useful tool for the recovery of this species, as well as for R. clavata, another elasmobranch with conservation concerns, currently classified as near threatened (IUCN, 2017). Positive responses to spatial management were previously reported for R. clavata in the Mediterranean by Dimech et al. (2008). Hunter et al. (2006) and Wiegand, Hunter & Dulvy (2011) also concluded that spatial closures were an appropriate approach for this species in the North Sea. Although the distribution of the local populations is probably wider than the MPA, many individuals might be using the area on a regular basis, given the known tendency to repeatedly return to specific locations (Hunter et al., 2005; Ellis et al., 2011).

Concerning C. lucerna, few studies have addressed its abundance in protected areas, despite its commercial importance. Piet & Rijnsdorp (1998) reported an increase in abundance of this species after the establishment of the “plaice box” in the North Sea, while Guidetti & Claudet (2010) found that it occurred in lower abundances inside the Torre Guaceto MPA (Italy) in comparison to the outside area. Our results show that this species can benefit from spatial closures. Within the Soleidae species (represented by nine species, four of which with frequency of occurrence >15%), M. azevia was the one presenting more positive indicators. It showed increases in abundance, biomass and size, and these trends may be related to the implemented regulations. As for the flounder C. linguatula, no trends over time were found in the study by Iannibelli & Musmarra (2008) in an area with trawling restrictions, while in our results, this was the only non-target species with consistent positive indicators. Claudet et al. (2010) pointed out that bycatch species may also be affected by fishing pressure or habitat degradation, which seems to be the case for this species.

The only species that exhibited a decrease in abundance was S. senegalensis, even though it was more abundant in the protected area and showed a size increase. This abundance decrease might be related with migrations to nearby estuaries (Andrade, 1990; Vinagre, Costa & Cabral, 2007) or to fishing grounds outside the reserve, with increasing exposure and vulnerability to gillnets and trammel nets. This phenomenon was hypothesized by Abecasis, Afonso & Erzini (2014) for this species in this same MPA (acoustic telemetry data).

It is also noteworthy that some species showed affinity to the BA, particularly P. lascaris and A imperialis. Both species had low frequency of occurrence and they exhibited restricted distributions in the study area. Similarly, C. linguatula was also unevenly distributed, but it showed affinity with the muddy bottoms in the FPA. This is in accordance with previous studies, that identified species with pronounced site-specific responses within MPAs (Claudet, García-Charton & Lenfant, 2011; Eddy, Pande & Gardner, 2014). One possibility is that their patchy distributions are related to biotic factors such as niche breadth, perhaps habitat preferences and/or prey availability. Studies such as the conducted by Ross (1986) have shown that among fish assemblages, the availability of food items is commonly an important factor for habitat discrimination. Moreover, it is argued that specialist species are more vulnerable to environmental stressors and exploitation, and that their dependence on habitat heterogeneity is also higher (Wilson et al., 2008; Slatyer, Hirst & Sexton, 2013).

As with rocky reef fish assemblages (Micheli et al., 2004; Ashworth & Ormond, 2005; Blyth-Skyrme et al., 2006), soft-bottom fishes exhibit a wide diversity of responses to spatial regulations. Our results suggest that part of this variability is connected with habitat heterogeneity. Still, the cause of site-specific preferences of some of the studied species is not fully understood. Studies of species movement ecology would provide insights for both habitat preferences and migratory behavior. Our data illustrates this gap of knowledge, as some mismatches related with migratory species were found. For instance, both R. clavata and R. alba demonstrate that species with yearly migrations can benefit from spatial measures, in accordance with the results of other studies (e.g., Claudet et al., 2010). However, it would be simplistic to infer that home range does not play a role in the effectiveness of spatial protection. Instead, possible connections are probably more complex, as illustrated by the species M. mustelus, a species that shares life history traits with the rays. This species did not exhibit consistent positive response to protection, and its migratory habits might contribute for this (strictly migrant vs seasonally migrant). The complexity of the effect of movement patterns in species response to protection was pointed out in previous studies (Claudet et al., 2010; Villegas-Ríos, Moland & Olsen, 2017).

In relation to species traits and their role in the response to protection, our most compelling results point out to body size as a possible factor. In our study, this is supported by the two rays (R. clavata and R. alba), two elasmobranchs with large body size and late maturity. This is in accordance with the observations of other authors (Claudet et al., 2006; García-Charton et al., 2008; Claudet et al., 2010; Sciberras et al., 2015), who reported body size as an important trait related to spatial protection. This might be due to the fact that larger fish may be more easily caught by fishing gear, and also because large body size generally reflects specific life history traits such as late maturity (Hutchings et al., 2012).

Implications and suggestions for future research

In this study, there was some evidence of protection effectiveness, supported by the significant increase in biomass in the protected areas. However, like in rocky habitats, effectiveness is cross-linked with habitat and species characteristics. In relation to habitat, it is important to identify habitat requirements of species that are protection targets, and habitat mapping should be carried out before the management plan development. This was also emphasized by other authors for rocky bottoms (e.g., Gaines et al., 2010; Halpern, Lester & McLeod, 2010). In addition, we recommend that habitat continuity should be taken into account in the design of MPAs aiming to protect soft-bottom fishes. Furthermore, we highlight the importance of enforcement measures, as compliance is essential for effective protection, long-term benefits and appropriate monitoring (e.g., Bergseth, Russ & Cinner, 2015). With respect to species characteristics, we highlight body size as a possible factor of response to protection. The inclusion of longer time scales would provide information on this and other underlying factors, and also reveal if more species could respond positively to spatial regulations. Future research on fish movements (telemetry), habitat requirements and trophic ecology (stable isotopes) would also offer useful knowledge to better understand how spatial protection measures should be used in the management and conservation of soft-bottom fish species.