The Influence of Boundary Habitat Continuity on Spillover from a Mediterranean Marine Protected Area

The effectiveness of marine protected areas (MPAs) to restore populations of exploited species, both within and outside of their boundaries through net movement of individuals (“spillover”), can potentially be affected by continuity of habitats across the boundaries. Sandy seabeds may reduce movement of reef-associated species across MPA boundaries, thereby increasing the ‘reserve effect’ while decreasing spillover. Underwater visual censuses were undertaken inside the Cerbère-Banyuls Marine Reserve (CBMR) (France) and adjacent non-protected areas to assess the influence of habitat on spillover. Total fish biomass and mean fish size were significantly higher within the MPA, but rapidly declined across the reserve boundary. Nevertheless, there was no indication of a sharper decline in biomass at the northern boundary where a habitat discontinuity was present relative to the southern boundary with continuous habitat. This result may reflect a number of complicating factors that make assessment of spillover potential difficult, and which may also lead to the uncertainty about which situations and how much spillover may contribute to fished populations outside reserves. In particular, the home range area of the key exploited species relative to the scale of the habitat mosaic, and potentially different levels of fishing pressure at each boundary likely contribute to variability. While the CBMR appeared particularly well-suited to investigating this question, resolving these issues and identifying general principles for where and how much spillover occurs will likely be difficult without a series of specially designed MPAs. This highlights a conundrum facing MPA establishment in the face of pressures to be successful for both biodiversity conservation and to offer fisheries benefits—the latter are clearly not ubiquitous, but a shortage of suitable MPAs that can be used as scientific tools for better understanding how and when these benefits may occur is precluded by a general lack of MPAs designed and managed for this purpose. The results of this study do, however, clearly highlight the biodiversity conservation benefits of the CBMR.


Introduction
The effects of marine protected areas (MPAs) on fish communities have been widely studied (Babcock et al. 1999;Planes et al. 2008), especially on species targeted by fishers . After the establishment of a MPA in an intensely-fished area, the biomass and body size of the target species generally increases through time within MPA boundaries (Halpern and Warner 2002;Harmelin-Vivien et al. 2008). Other metrics like fish abundance and species richness may also increase (Vanderklift et al. 2013), although such effects appear less general and more idiosyncratic than recovery of exploited species biomass, and are generally not considered as unequivocal indicators of MPA effectiveness Lenihan et al. 2021).
MPAs exclusively protect those species that remain inside the protected area (Chapman and Kramer 2000). Many temperate rocky reef species are site-attached and show some degree of fidelity to certain habitats (García-Charton and Pérez-Ruzafa 2001). However, random movements, relocation of home ranges (Kramer and Chapman 1999) or 1 3 increasing densities of large individuals within the MPA can result in emigration of adults to adjacent fished locations, a process contributing to spillover (along with juvenile and larval dispersal outside of boundaries; Stobart et al. 2009;Di Lorenzo et al. 2016).
Spillover from MPAs may occur in relatively mobile or large species which are less site-attached, or as a result of density-dependent or competitive interactions within highly effective MPAs. The type of habitat bordering the MPA influences movement of species, so can potentially influence the degree of spillover from MPAs (Barrett et al. 2007;Forcada et al. 2008). Continuous habitats should facilitate spillover (Goñi et al. 2008), while habitat discontinuities may form natural boundaries and reduce emigration rates for species with particular habitat preferences (Barrett 1995;Kramer and Chapman 1999;Chapman and Kramer 2000).
The responses of fish assemblages to protection in MPAs can be difficult to differentiate from the effects of a number of environmental and habitat-related factors (Jaworski et al. 2010). Some factors, such as habitat structure, have been demonstrated to have stronger effects on fish assemblage structure than fishing regulations (García-Charton et al. 2004). Depth and recruitment potential have also been shown to influence potential recovery of fish populations within MPAs, and detectability of such effects (Claudet et al. 2006;Crec'hriou et al. 2008;Pastor et al. 2009;Valle and Bayle-Sempere 2009).
The aim of this study was to assess spillover patterns in fish biomass from a Mediterranean MPA to evaluate the existence of a reserve effect, and test the hypothesis that fish spillover is lower across a boundary characterized by discontinuous habitat than a boundary characterized by continuous habitat. The study capitalizes on the distribution of rocky reef habitat along the south-western corner of France, where the Cerbère-Banyuls Marine Reserve (CBMR) provides a natural experiment involving removal of fishing pressure, punctuated by contrasting habitat continuity along a linear coastline. The habitat connection at the northern boundary of the CBMR is disrupted by a 400 m-long sand flat, while the southern boundary presents a continuous fringing rocky reef. Considerable water depth offshore (> 60 m) likely limits adult movements of most shallow water species to the linear coastal reef habitat.
The three key elements to the study were: (i) to test whether reserve protection has benefited the fish community; (ii) to test for differences in fish biomass between fished locations positioned across the northern and southern boundaries, which may relate to different spillover patterns due to differences in habitat continuity at the boundaries of the reserve, (iii) to explore the importance of habitat structure and depth in contributing to observed patterns.

Site Description
The Cerbère-Banyuls Natural Marine Reserve (CBMR) extends along 7 km of shoreline and out to 1.5 nautical miles from shore (Pastor et al. 2009) in the Northwest Mediterranean region, at the French province of Rousillion. The CBMR was legally created in 1974 (Pastor et al. 2009;Planes et al. 2008), and in 1979 an Integral Reserve area was established within the Protected zone. The CBMR covers an area of 650 ha separated in two zones with different levels of protection ): (i) Integral Reserve Zone, which covers an area of 65 ha, in which professional fishing, SCUBA diving, angling, spear fishing and anchoring are forbidden; (ii) Protected Zone, which covers ca. 585 ha and in which spear fishing is forbidden, and fishing, anchoring and cruising are regulated.
The CBMR area is dominated by rocky substrates, with medium to large boulders (> 2 m diameter), vertical walls and gullies located predominantly inside the reserve. However, at the boundaries of the protected area the habitat distribution differs along the Reserve margin; to the south, fringing rocky reef extends relatively continuously across the Reserve boundary ( Fig. 1), while to the north, a sandy embayment and long sand beach extends from the boundary of the protected area, separating the rocky habitats inside the reserve from those located at the adjacent non-protected areas by a minimum of 400 m.
Sampling For the study, 13 sites were selected according to their protection status (Fished: stations S1, S2, S3, S6, S8, S9, S10 and S11; Boundary (fished-Protected): S7 and S13; Protected: S5 and S12, and Integral: S4) ( Fig. 1) and distance from the centre of the reserve, with a greater survey intensity across MPA boundaries in both directions along the coast. At each site, two 50 m length transects were laid from a given point in opposite directions parallel to the shore over rocky substrate and at constant depth.
At each site, the fish community was assessed using underwater visual census along 50 m transect lines × 10 m wide (total area surveyed 500 m 2 per site), following Reef Life Survey methodology (Edgar and Stuart-Smith 2014; and described in detail in an online methods manual at https:// reefl ifesu rvey. com). Two surveys were consistently undertaken at each site, starting from a given point and running in opposite directions parallel to the shoreline in a depth range from 4 to 11 m (Table 1) in order to minimize the influence ◂ of depth on fish community structure. The average depth of transect lines was recorded, and used as a covariate in analyses. Likewise, habitat complexity was estimated for each transect using a semi-quantitative technique (Stuart-Smith et al. 2008). This consisted of two scores allocated for each transect to account for horizontal and vertical structure of the seabed. Horizontal structure was defined by the presence of caves, holes or crevices present along the 50 × 10 m survey area, using three categories: (i) numerous caves and overhangs (> 50% coverage); (ii) some caves or overhangs (10-50% coverage); (iii) no caves or overhangs (< 10% coverage). Vertical relief was estimated as the relative depth variation from the shallowest to the deepest point of the transect, together with the presence of large scale structures such as walls and gullies. The three categories used to score vertical structure were: (i) high relief (presence of boulders greater than 5 m high, walls or gullies; high variation in profile along the transect); (ii) medium relief (presence of small walls and/ or few small boulders; moderate profile variation along the transect); (iii) flat (flat or gently sloping area, without any large structures or variations in profile along the transect).
Distance from the Reserve (DCCentre) was calculated using ArcGIS software. This was the closest distance by water between each survey site and the geographical "centre" of the Reserve (defined as the midpoint between the boundaries of the Integral zone), measured following the depth contour and avoiding crossing sand extensions. This effectively simulated the movement of a fish around the coast (Bell 1983) when following rocky reef habitat (Chapman and Kramer 2000;Forcada et al. 2009).
Fish biomass was calculated using the length-weight relationship: W = aL b (Froese and Pauly 2000), where W is derived weight in grams, L is fish Total Length (TL) in cm, and parameters a and b obtained for each species from published relationships for Mediterranean species (Valle et al. 2003;Froese and Pauly 2011). When these parameters were only available for relationships using other length measures (standard length or fork length) they were converted to the appropriate length measure using published length-length relationships (Froese and Pauly 2011). Size estimates made by divers were first transformed following published size corrections (Edgar et al. 2004).

Statistical Analysis
Generalised linear models (GLMs), using a negative binomial family-error structure because of overdispersion resulting from a large number of zeros, were fitted to both matrices of fish abundance and fish biomass via the R package 'mvabund' (Wang et al. 2012). GLMs were used to analyze fish biomass (total biomass of fishes per 500 m2) and fish size (mean size of individuals per 500 m2) associated with the protection status of sites: (i) Fished: sites located outside the reserve, (ii) Boundary Protected-Fished: sites located at the boundaries between the protected zone and the fishing zone, (iii) Protected: sites located inside the protected zone and, (iv) Integral: sites located at the boundary of the integral reserve zone. The second was used to test for significance of differences in total fish biomass between (i) fished sites located outside the continuous habitat (southern) boundary of the reserve (S6, S9, S10, S11), vs. (ii) fished sites located outside the discontinuous habitat (northern) boundary of the reserve (S1, S2, S3, S8). A distance-based redundancy analysis (db-RDA, Legendre and Anderson 1999) tested whether variation in any of the environmental variables significantly contributed to explain variation in the fish community structure in areas with different degrees of protection. Multivariate multiple regression, using the DISTLM (DISTance based Linear Model) routine (Anderson 2001) was also used to assess the contribution of the distance from the centre of the Reserve on total fish biomass, after accounting for the natural variation resulting from depth and habitat complexity. The AIC routine was used as a selection criterion (Legendre and Anderson 1999) and depth and habitat complexity were forced inclusions in the DISTLM, with distance (DCC) added to test for significance and unique contribution to the variance in fish biomass over and above the effects of depth and habitat complexity. Lastly, non-parametric multidimensional scaling (n-MDS) was undertaken on the multivariate fish community structure data, with biomass values for each species at each site logtransformed and used to generate a Bray-Curtis Dissimilarity matrix on which the n-MDS was performed. The DistLM and n-MDS were carried out using the PRIMER v6 and PER-MANOVA + statistical package.

Results
A total of 13,000 m 2 of seabed was surveyed across 13 sites, and data from 3,184 individuals of 24 fish species recorded. Total fish biomass averaged 70.2 kg per 500 m 2 .

Fish Biomass and Size
A clear distinction was evident in total fish biomass between sites in the different protection classes (Fig. 2). Highest biomass was observed at the integral site (21.3 kg /500 m 2 ), followed by protected sites (10.3 ± 0.4 kg /500 m 2 ), then sites located at the boundary between protected and fishing zones (6.9 ± 0.5 kg /500 m 2 ), while fished Fig. 2 Total biomass of fishes (kg /500 m 2 ) for fished sites (S1, S2, S3, S6, S8, S9, S10, S11), sites located at the external protected area boundary (S7, S13), protected sites (S5, S12) and one site located at the boundary of the integral reserve zone (S4). Error bars show standard error of the mean sites showed the lowest biomass values (1.8 ± 0.2 kg /500 m 2 ). However, this difference in fish biomass was similar between fished sites located at the continuous boundary (1.7 ± 0.5 kg /500 m 2 ) and those located at the discontinuous boundary (1.9 ± 0.1 kg /500 m 2 ) (Fig. 3).
Mean fish size was greatest at the integral site (17.9 cm), followed by the protected sites (15.9 ± 0. 8 cm), external boundary sites (15.4 ± 0.1 cm), and fished sites (12.1 ± 0.4 cm), which hosted the lowest values (Fig. 4). Within the fished sites, those located outside the continuous boundary showed slightly lower mean fish size (11.9 ± 0.8 cm) than those located outside the discontinuous boundary (12.3 ± 0.5 cm), but these differences were not significant (GLM, estimate = 0.1, p = 0.29). Significant differences in mean size were found between fished sites Fig. 3 Total biomass of fishes (kg /500 m 2 ) for fished sites outside the discontinuous boundary (S1, S2, S3, S8) and the continuous boundary (S6, S9, S10, S11). Error bars show standard error of the mean

Fishing
Boundary After accounting for the expected contributions of depth and habitat complexity on total fish biomass (Table 2), the distance following the depth contour to the centre of the Reserve (DCCentre) added a significant contribution to the variance in fish biomass explained (11%; Table 2; Fig. 5). The best model therefore included DCCentre and explained 43% of total variation in mean fish biomass.

Influence of Protection on Fish Assemblages
Protection effects were most noticeable in elevated biomass of several exploited species, i.e. Diplodus cervinus, D. puntazzo, D. sargus and D. vulgaris, and Sarpa salpa. The n-MDS showed clear patterns in fish community structure among all site groups (fished, external boundary, protected and integral Reserve) (Fig. 6). The similarity between fished and protected sites was 52%, and within each of the groups (fished and protected) the similarity of fish community structure was 60%. Table 2 DISTLM results for testing the effects of distances of sites from each boundary on total fish biomass. Depth and habitat complexity were forced inclusions, together explaining 29.0% of variation in fish biomass. DCCentre = distance following the depth contour to the centre of the reserve. Bold letters denote significant differences. The final model, including Depth, Habitat complexity and DCCentre explained 43% of variation in mean fish biomass

Fig. 5
Distance-based redundancy analysis (db-RDA) biplot of first and second axes for fish biomass at sites located within and outside the MPA. A vector overlay indicates relationships with covariates: DEPTH, average depth at the sampled site; DCCentre, distance following the depth contour to the centre of the reserve; HABCOM, habitat complexity

Discussion
We found no clear support for our prediction that discontinuous reef habitat at the reserve boundary would limit fish dispersal, and be reflected in a more pronounced gradient in fish biomass across the boundary. To the contrary, fish biomass did not differ significantly between the northern fished sites located outside the discontinuous boundary and southern sites connected to the reserve by a continuous fringing reef. Significant differences in total fish biomass were evident between sites at the boundary of the protected zone, sites inside the protected zone, and those at the boundary of the integral reserve. This outcome is consistent with findings of previous studies, where higher biomass was present within the protected areas of Mediterranean MPAs than in the surrounding non-protected zones (Harmelin-Vivien et al. 2008;García-Charton et al. 2008). The order of magnitude greater biomass of fish at the integral site inside the CBMR in this study compared to adjacent fished sites strongly suggests that either this MPA is highly productive and can maintain high fish biomass in the face of considerable spillover, or that relatively limited dispersal of adults occurs across the reserve boundaries, regardless of habitat continuity. The dominant large fishes in this MPA (and most Mediterranean MPAs) are sparids, which are generally known to have high site fidelity and to respond well to protection in MPAs (Kerwath et al. 2007(Kerwath et al. , 2013. Despite being a small MPA, CBMR is relatively old by global standards, hence the fish community within the reserve has had a long time to recover from the very high fishing pressure experienced prior to protection, and that is still occurring outside the reserve. Our results are not likely to be confounded by differences in habitat complexity or depth between the reserve and fished sites, two factors that can affect the structure of fish assemblages and explain small-scale spatial patchiness in fish communities (García Charton and Pérez-Ruzafa 2001;Harmelin-Vivien et al. 2008). Our DistLM on fish biomass data included depth and habitat complexity as forced inclusions, thereby specifically testing whether the distance of sites from the centre of the reserve helped explain the remaining variation in fish biomass between sites. The decline in biomass from the integral reserve to the outside fished areas was still evident after accounting for variation in depth and habitat complexity, supporting that this trend was a strong effect of protection.
Analyses of fish size generated results consistent with those for overall biomass, providing further support to Fig. 6 n-MDS representation of fish community structure in relation to its protection class (fished, external boundary, protected, integral). Square root transformed biomass data for each species were used and Bray-Curtis dissimilarity 1 3 significant protection effects, albeit with similar values at sites along the continuous and discontinuous habitat boundaries. Thus, the elevated biomass within the reserve is at least in part driven by an increase in body size, regardless of abundance changes, as has been found in previous MPA research (Babcock et al. 1999).
The CBMR appeared particularly well-suited to investigating the influences of habitat discontinuities on spillover at a scale that is feasible to assess using visual censuses, and in a situation in which spillover has been previously confirmed from studies of fisheries effort and catches (Goñi et al. 2008(Goñi et al. , 2010. While no effect of habitat continuity was found in our study, results of other studies suggest variability in such an effect (Di Lorenzo et al. 2016;Forcada et al. 2008Forcada et al. , 2009). This variability likely arises from numerous factors such as each species' home range size and habitat specificity, the scale of habitat patchiness and fishing pressure and access at each boundary. Indeed, a habitat effect is likely to vary as much as the presence of a spillover effect in general. Understanding in which situations it is likely important, and what trade-offs exist between biodiversity conservation outcomes inside the reserve vs fisheries production outside of the reserve will require much more study, and possibly experimental designs that are not presently possible. Moving forward in this field would benefit greatly from a globally coordinated series of MPAs designed and managed as scientific tools for understanding these trade-offs and general principles. In the absence of such a network, addressing short-comings in management effectiveness of a majority of the world's MPAs (e.g. due to compliance, no-take regulations, size, age; Gill et al. 2017) will assist in boosting the value of these areas for undertaking the science needed to support the design and roll-out of more effective MPAs in future.
The present study provides only a preliminary investigation of the specific role of habitat continuity on fish spillover to nearby fished sites outside protected areas. An integrative study, incorporating other methodologies and techniques such as tagging, environmental DNA (eDNA), landings of exploited species and fishing effort data is needed to provide an improved picture on the functioning of protected areas and nature of spillover effects into adjacent zones. Importantly, a more general understanding of the processes that lead to and promote spill-over will require coverage of many reserves with varying levels of habitat continuity, different fish community compositions and environmental settings.
Acknowledgements This study was supported by the European Community-ASSEMBLE grant agreement nº. 227799. We are grateful to Dr. Hilmar Hinz, Dr. John Turner and Dr. Andrew Davies (University of Bangor, Wales) for their advice, support with GIS methodology and diving protocols. To Dr. Jemina Stuart-Smith for her assistance throughout the field surveys. Finally, to the staff of Laboratoire Arago-Observatoire Océanologique de Banyuls (France) and especially to Dr. Pascal Romans for the assistance provided with the surveys, sites information and species guidance. To Dr. R. Stuart-Smith and G. Edgar (University of Tasmania) for their support and encouragement throughout these years and also for improving an earlier version of this manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Statement of Conflict of Interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.