Mediterranean marine protected areas have higher biodiversity via increased evenness, not abundance

1School of Zoology, George S. Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv, Israel; 2German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany; 3Department of Computer Science, Martin Luther University Halle-Wittenberg, Halle (Salle), Germany; 4Stazione Zoologica Anton Dohrn, Dipartimento Ecologia Marina Integrata, Sede Interdipartimentale della Sicilia, Lungomare Cristoforo Colombo (complesso Roosevelt), Palermo, Italy; 5Consorzio Interuniversitario per le Scienze del Mare, CoNISMa, Rome, Italy; 6Université Côte d’Azur, CNRS, UMR 7035 ECOSEAS, Nice, France; 7Department of Biology, University of Vermont, Burlington, VT, USA; 8Institute of Biology, Martin Luther University Halle-Wittenberg, Halle (Saale), Germany; 9Department of Community Ecology, Helmholtz Centre for Environmental ResearchUFZ, Halle (Saale), Germany; 10Leuphana Universität Lüneburg, Lüneburg, Germany; 11Biology Department, College of Charleston, Charleston, SC, USA; 12Hopkins Marine Station and Stanford Center for Ocean Solutions, Pacific Grove, CA, USA; 13National Geographic Society, Washington, DC, USA and 14The Steinhardt Museum of Natural History, Tel Aviv University, Tel Aviv, Israel


| INTRODUC TI ON
Protected areas are important for conservation strategies in marine and terrestrial systems (Gaston, Jackson, Cantú-Salazar, & Cruz-Piñón, 2008;Watson, Dudley, Segan, & Hockings, 2014). They protect biodiversity by reducing mortality due to habitat destruction and harvesting. Abundance and biomass are often higher inside protected areas (Coetzee, Gaston, & Chown, 2014;Edgar et al., 2014), whereas empirical evidence for biodiversity gains within protected areas is mixed (Gaston et al., 2008;Lester et al., 2009). Marine protected areas (MPAs) are often designed and implemented for a combination of biodiversity conservation and to support sustainable fisheries (Gaines, White, Carr, & Palumbi, 2010), and studies examining protection effects on biodiversity typically quantify species richness changes at the scale of individual protected areas (White et al., 2011). However, quantifying species richness at a single scale provides an incomplete picture of how biodiversity changes in response to an external driver (e.g. Chase & Knight, 2013;Hillebrand et al., 2017;Supp & Ernest, 2014). Moreover, MPAs impose spatial variation in exploitation, and are often part of protected area networks (Wood, Fish, Laughren, & Pauly, 2008). A multi-scale approach is needed to more fully evaluate the influence of protection on patterns of biodiversity.
Increased fish abundance and biomass are the strongest and most commonly observed responses to protection inside MPAs (Lester et al., 2009;Soykan & Lewis, 2015). Species richness is also often greater inside MPAs (Lester et al., 2009), although gains are typically smaller relative to those of biomass and abundance (Soykan & Lewis, 2015). Here, we examine the multiple pathways that influence how species richness increases with sampling effort and spatial scale (i.e. the species accumulation curve). The species accumulation curve is known to be influenced by three components of the underlying community: changes in the number of individuals, changes to the relative abundance of species and/or changes to patterns of spatial aggregation (Chase & Knight, 2013;He & Legendre, 2002;McGill, 2011). As a result, it is useful to explore how these underlying components change and contribute to biodiversity patterns across scales (Chase et al., 2018;McGlinn et al., 2019).
Protection from harvesting inside MPAs potentially affects all three components underlying species richness and its scaling (Tittensor, Micheli, Nyström, & Worm, 2007). Higher abundances of species targeted by fisheries are one of the most commonly observed responses to protection (Claudet et al., 2010), and communities with more individuals typically have more species via the 'more individuals hypothesis' (Storch, Bohdalková, & Okie, 2018).
Increased abundances of fishery target species-which usually occupy high trophic levels (Pauly, Christensen, Dalsgaard, Froese, & Torres, 1998)-may also alter the overall evenness of the community. Increasingly abundant predators may influence the total and relative abundance of prey, possibly reducing the overall variability among species abundances (Soykan & Lewison, 2015). Thus, quantifying changes to patterns of commonness and rarity among species is required to understand the response of biodiversity to protection.
Moreover, MPA networks introduce site-to-site variation in protection from exploitation. If this spatial variation in exploitation changes spatial patterns of within-species aggregation (Baskett & Barnett, 2015), then this too may alter biodiversity at local and regional scales (McGill, 2011).
In addition to introducing spatial heterogeneity in protection from exploitation, MPA networks are sometimes designed to maximize complementarity (i.e. the diversity accumulated across sites; Margules & Pressey, 2000). For example, some planned MPA networks accumulate diversity across sites by protecting different habitat types or by incorporating different human use regulations (e.g. Fernandes et al., 2005). This further emphasizes the need for assessments of MPA networks at both local and regional the regional scale than at the local scale, meaning that relative abundances at the regional scale were less even than at the local scale.

Synthesis and applications.
Although marine protected areas (MPAs) are known to strongly influence fish community abundance and biomass, we found that changes to the relative abundance of species (i.e. increased evenness) dominated the biodiversity response to protection. MPAs had more relatively common species, which in turn led to higher diversity for a given sampling effort. Moreover, higher β-diversity of common species meant that local-scale responses were magnified at the regional scale due to site-to-site variation inside protected areas for exploited species. Regional conservation efforts can be strengthened by examining how multiple components of biodiversity respond to protection across spatial scales.

K E Y W O R D S
beta-diversity, biodiversity, conservation, marine protected areas, protected areas, scale dependence spatial scales (Grorud-Colvert et al., 2014;Socolar, Gilroy, Kunin, & Edwards, 2016;White et al., 2011). β-diversity, the component of regional biodiversity (γ-diversity) that describes the between-site differences in the diversity of local assemblages (α-diversity), should provide information of how local responses combine at the network scale (Socolar et al., 2016). For example, if the same, formerly exploited species returned to all protected sites within an MPA network, site-to-site variation would likely decrease within reserves (additive homogenization; Socolar et al., 2016).
Alternatively, if different exploited species returned to different protected sites inside an MPA network, as might be expected if sites were selected to maximize habitat diversity, then β-diversity would be expected to increase (additive differentiation ;Socolar et al., 2016).
We examine how a regional system of MPAs in the Mediterranean Sea affects fish biodiversity and its scaling. Coastal regions of the Mediterranean are home to more than 150 million people, and multiple human stressors have impacted ecosystems for centuries (Guidetti et al., 2014;Micheli et al., 2013). Currently, 6.5% of the Mediterranean Sea is designated with some level of protection, and 0.04% is fully protected (PISCO & UNS, 2016). We evaluated how fish biodiversity across multiple scales responds to protection by dissecting species richness into components: the number of individuals, the relative abundance of species and the patterns of within-species aggregation. Examining the responses of multiple biodiversity components across scales reveals new insights into how fish communities respond to protection.

| Data
Fish assemblages were sampled in the northern Mediterranean Sea ( Figure 1) during May-June 2007 and 2008 (see Guidetti et al., 2014;Sala et al., 2012 for further details). At each site, similar depths (8-12 m) and habitats (rocky reefs) were selected to minimize environmental heterogeneity, and fish assemblages were surveyed using three 25 × 5 m strip transects.
Our regional-scale analyses required a sampling design where fished and protected sites encompassed similar spatial extents (e.g. to minimize influences other than protection on β-diversity).
Accordingly, we reduced the extent from Guidetti et al. (2014), and grouped MPAs classified as having high or intermediate protection as 'protected' (n = 43 protected sites, representing eight marine protected areas), and non-enforced MPAs and fished sites as 'fished' sites (n = 41 fished sites, representing seven fished areas).
We examined the sensitivity of our discrete-scale analyses to the simplified protection classifications by separating 'protected' sites into fully-(well enforced, no-take MPAs; n = 21) and partially protected sites (i.e. MPAs where some fishing is allowed or some illegal fishing may occur due to weak enforcement; n = 22; Figure S1). However, because our multi-scale analyses rely on pairwise comparisons between rarefaction curves, we chose to use only protected and fished categories to simplify the presentation of results.
Our samples from fished and protected sites were not matched spatially (i.e. we do not have samples from inside and outside protected areas at all locations; Figure 1). To make inferences of protection effects as robust as possible, we gathered additional data to adjust for variation in the environment and other human impacts.
Habitat complexity was measured in situ along each transect as substrate rugosity (see Guidetti et al., 2014), and environmental covariates (e.g. temperature, chlorophyll A concentration, etc.) were extracted from Bio-ORACLE (Tyberghein et al., 2011, see appendix A for details). For a proxy of human pressure, we used the cumulative impact layer that integrates 22 anthropogenic drivers (e.g. various types of fishing, invasive species, climate change, nutrient input) for the Mediterranean from Micheli et al. (2013). These data were included as covariates in our α-scale analyses, and we used permutation tests to examine for systematic differences between fished and protected areas across all sites. All covariates were mean centred and standardized by dividing by one standard deviation prior to all analyses. PERMANOVA on a Euclidean distance matrix did not reveal strong evidence for systematic differences between fished and protected sites (F = 2.2, p = .08; Table S1), and the variance of the covariates did not differ between fished and protected sites (F = 2.5, p = .14; Table S2).
To assess whether our results were likely to be strongly influenced by missing (i.e. unobserved) species, we calculated abundance-based coverage (Chao & Jost, 2012 Figure S2), meaning the probability that another individual sampled would represent a new species was <2.5%.

| Biodiversity dissection and scale dependence
We examined the scale-dependent response to protection using complementary discrete-and multi-scale analyses (
We calculated species richness, total number of individuals and a measure of species relative abundances at the α-and γ-scales.
Examining the total number of individuals (N) provides insight into whether richness changes are simply due to different numbers of individuals being sampled. To assess whether changes in relative abundance were underpinning altered species richness, we calculated the probability of interspecific encounter (PIE). The PIE is the probability that two individuals sampled randomly from a community are of different species (Hurlbert, 1971), and higher values represent more even communities. We transformed the PIE into an effective number of species (S PIE ) that has the same units as species richness (Jost, 2006). Finally, we calculated species richness (S) and rarefied species richness (S n ; expected richness for n individuals, Gotelli & Colwell, 2001). S is more sensitive to rare

Metric Definition Interpretation
Discrete scale N Total number of individuals Measure of how the density of individuals responds to protection. N scales approximately linearly with area (i.e. N is scale independent) and so we only calculated N at the local (α) scale α S PIE , γ S PIE Number of equally abundant species needed to yield the observed Probability of Interspecific Encounter (PIE, Jost, 2006). Equivalent to diversity of q = 2 (Jost, 2007) Differences in S PIE reflect changes in the effective number of relatively common species (Jost, 2007), or equivalently due to the relationship with the PIE, changes in evenness α S n , γ S n Expected number of species for n individuals (Hurlbert, 1971); calculated at the α-and γ-scales Differences in S n reflect changes in the SAD a only, the effects of aggregation and N are removed α S, γ S Observed species richness at the scale of sites ( α S), or all fished or protected sites combined ( γ S) Differences in S are due to some combination of changes in N, the SAD and/or within-species aggregation β-S PIE Ratio of γ S PIE over average α S PIE Number of distinct communities at the regional scale. Higher values of β-S PIE reflect greater site-to-site variation mostly due to aggregation of common species β-S Ratio of γ S over average α S Number of distinct communities at the regional scale. Higher values of β-S reflect greater siteto-site variation due to changes in N, the SAD and aggregation of common and rare species Multi-scale SAD effect Calculated as the difference between the individual-based rarefaction curves Quantifies the contribution of changes in the SAD to observed changes in species richness continuously across scales N effect Calculated by subtracting the difference between the individual-based rarefaction curves (SAD effect only) from the difference between the two non-spatial curves (N and Quantifies the contribution of changes to patterns of within-species aggregation to observed changes in species richness continuously across scales species, whereas S PIE depends on the number of abundant (or common) species (Jost, 2006), and combined they provide complementary information on how rare and common species respond to protection. Additionally, comparisons of changes in S with changes in S n reveal whether changes in the number of individuals (N) are contributing to diversity patterns. For example, if protection effects on S are not found on S n , then changes in N dominate the gains in species richness. However, if protection effects on S and S n are found, then changes in both N and the species abundance distribution (SAD) are contributing to the biodiversity response (Chase et al., 2018).
Protection effects on biodiversity at the α-scale were quantified using hierarchical linear models. The total number of individuals (N), S PIE and S n were modelled assuming log-normal distributions and an identity-link function; species richness (S) was modelled assuming a Poisson distribution and a log-link function. All models included the environmental and cumulative human impact as continuous (mean centred and standardized) covariates; status as fished or protected was coded as a categorical covariate; sites were grouped into the distinct protected and fished areas that they came from, and this location was included as a random intercept. For Bayesian inference and estimates of uncertainty, models were fit using the Hamiltonian Monte Carlo (HMC) sampler Stan (Carpenter et al., 2017), and coded using the 'brms' package (Bürkner, 2017). All models were fit with four chains and 2,000 iterations, with 1,000 used as warmup. We used weakly regularizing priors and visual inspection of the HMC chains showed excellent convergence.
At the γ-scale, our design was unbalanced (41 fished and 43 protected sites), so comparisons were made with sample-based rarefaction (Gotelli & Colwell, 2001). Using bootstrap resamples without replacement, 35 random sites were sampled from both fished and protected area 200 times. Species counts were accumulated within each treatment, and we calculated S PIE and S for each of the bootstrap resamples. The effects of protection at the γ-scale were examined using the median and the 95% quantiles of the resamples.
We also used these same bootstrap resamples to examine protection effects on β-diversity. For each of the resamples, we calculated α-scale means of S PIE and S. β-diversity (i.e. β-S PIE and β-S) was then calculated as the ratio of the resampled γ/α metrics. Similar to our γ-scale results, we examined protection effects using the median and 95% quantiles of all the resampled β-diversities.

| Multi-scale analyses
We

| Sensitivity to exploitation and the effects of protection
To examine whether the biodiversity response to protection depends on species sensitivity to exploitation, we retrieved a 'sensitivity to exploitation' score from FishBase (called 'vulnerability' in FishBase, but referred to hereafter as sensitivity; Froese & Pauly, 2017). Sensitivity is a continuous variable between 0 and 100, calculated using eight life-history traits (Cheung, Pitcher, & Pauly, 2005), where high scores represent high sensitivity to exploitation.
We performed all analyses for the whole community combined (total species richness, S = 51), and separately for fishes with high and low sensitivity to exploitation. We defined high and low sensitivity as the upper 30% (S high sensitivity = 16) and lower 70% (S low sensitivity = 35) quantiles of the sensitivity scores respectively; and examined whether our results remained qualitatively consistent when different quantile thresholds were used to define high and low sensitivity.
All data manipulation and analyses used r (R Development Core Team, 2017). FishBase trait data were accessed using rfishbase (Boettiger, Lang, & Wainwright, 2012), and the multi-scale analyses   Figure S3). When the protected sites were subdivided into those with partial and full protection ( Figure S1), results were qualitatively consistent: fully protected sites had more even relative abundances and more species than partially protected sites ( Figure S4).

| Discrete-scale analyses
At the γ-scale, the number of common species (S PIE ) increased with protection irrespective of sensitivity to exploitation (Figure 3a,c,e). Protection also increased the species richness for all fishes combined (Figure 3b) and exploited species (Figure 3d), but not fishes less sensitive to fishing (Figure 3f). These results were largely qualitatively consistent when the protected sites were divided into full and partial protection: the number of common species (i.e. evenness) increasing from fished through partially-to fully protected sites ( Figure S5a,c,e); however, partially and fully protected areas had similar species richness at the γ-scale ( Figure S5b,d,f).
Protection increased the β-diversity of relatively common species (β-S PIE ; Figure 4a,c,e), but there was no effect on β-S (Figure 4b,d,f).
β-S PIE was >1 for fishes most sensitive to exploitation, suggesting that there was more than one distinct community of these fishes at the regional scale. In contrast, β-S PIE values were <1 for low sensitivity species and all fishes combined. This means that there were fewer common species at the regional scale than the average local site, or equivalently, evenness was lower at the regional-compared to the local scale. Dividing protected sites into full and partial protection revealed β-S PIE was lowest among unprotected sites and approximately equal for partially and fully protected sites, and that β-S was highest among partially protected sites ( Figure S6).

| Multi-scale analyses
For all fishes combined, changes in the SAD made the largest contribution to species richness gains inside protected areas ( Figure 5). For fishes highly sensitive to exploitation, N and the SAD made scale-dependent contributions to richness gains: increased numbers of individuals (N) contributed most at intermediate scales ( Figure 5b), whereas the SAD contribution was largest at the full extent of the study (Figure 5b). This suggests that protection is influencing the whole species abundance distribution of highly sensitive species: there are more relatively common species in protected areas at local scales (Figure 2d), and rare species were accumulated across the extent of MPAs network (see absence of asymptote for the SAD effect in Figure 5b). Finally, the SAD contribution to gains of species least sensitive to exploitation was also scale-dependent, peaking at intermediate scales and returning to zero at the extent of the study ( Figure 5c). These results remained qualitatively consistent when we varied the cut-off value used to determine fishes with high and low sensitivity to exploitation ( Figure S7).

| D ISCUSS I ON
It is well-established that MPAs have positive effects on marine ecosystems (e.g. Edgar et al., 2014;Mellin, Aaron MacNeil, Cheal, Emslie, & Julian Caley, 2016). But pathways through which MPAs protect local and regional biodiversity are less well known. The simplest and most intuitive effect would be that increased fish abundances due to protection from harvesting lead to increased species richness in MPAs via the 'more individuals hypothesis'. That is, with more individuals in protected areas, we would expect more species via random sampling alone. While this effect certainly plays a role, we found that variation in numbers of individuals actually contributed very little to changes in species richness under protection, because numbers of individuals of all species combined did not vary between areas of different protection status ( Figure S3).
Instead of species richness changes due to altered numbers of individuals, we found that rare and common species were disproportionately affected by protection. Specifically, increased numbers of common species (or equivalently, increased evenness) was the most consistent biodiversity response inside protected areas at both the local (α) and regional (γ) scales, and our continuous analysis showed richness gains in protected areas for the whole fish community were largely due to changes in the SAD.
At the local (α) scale, these results are consistent with a recent

(i)
Protection meta-analysis of community-level MPA effects showing increased evenness in species' relative abundances (Soykan & Lewison, 2015). Here, we additionally show how altered patterns of relative abundance can make scale-dependent contributions to biodiversity gains.
Separating species into groups more or less sensitive to exploitation revealed distinct patterns. Exploited species responded most strongly to protection; richness gains at smaller scales were due to a combination of increased numbers of individuals and evenness, but the individuals' (N) effect was tempered, while the evenness (SAD) effect increased, with increasing scale. In contrast, species less sensitive to exploitation had smaller richness gains at small scales only, which were due solely to increased evenness.
The sensitivity to exploitation metric we used combines life-history traits known to influence vulnerability to exploitation (e.g. maximum size, age at first maturity; Cheung et al., 2005), and has proved successful at predicting population status without formal stock assessments (Reynolds, Dulvy, Goodwin, & Hutchings, 2005). However, where site or geographical variation in fish traits exists (Claudet et al., 2010), or where variation in MPA size, shape, fishing effort or gears (e.g. due to local regulations) alter species' exposure to fishing, local knowledge may be needed to better determine siteand species-specific vulnerability to exploitation (e.g. Claudet et al., 2010  Protection an assessment of European MPAs that used expert opinion to derive species-and location-specific sensitivities (Claudet et al., 2010).
What mechanisms could underlie changes in species' relative abundances following protection? Similar to existing work describing reduced total abundance of prey species within MPAs (Cheng, Altieri, Torchin, & Ruiz, 2019;Claudet et al., 2010), we hypothesize that increased evenness among prey species could reflect stronger top-down control inside protected areas. We found a trend towards lower abundances of prey species (see Figure S8 for the positive relationship between sensitivity and trophic level) in MPAs, and lower abundances were accompanied by increased evenness (and a reduction in the relative abundance of the most common prey species of almost 10%; Figure S9).
These results are consistent with density-dependent immigration out of MPAs; or predators focusing on the most common prey species, whereby increasing predator abundance following protection may disproportionately affect densities of abundant prey species.
Changes in evenness inside protected areas showed important scale dependence and spatial variation associated with sensitivity to exploitation. The increase in the number of common species with high sensitivity to exploitation was greater at the regional compared to local scales, reflecting a tendency for different protected sites to have different species. In contrast, fewer common species less sensitive to exploitation were found at the regional scale than at the typical local site, meaning the regional community had a less even SAD than the average local site for these species. The finding of increased β-S PIE for exploited species inside protected areas

(f)
Protection suggests protection could act to reverse taxonomic homogenization possibly associated with harvesting, and shows that local conservation initiatives can combine synergistically across a regional system of MPAs.
Not all marine protected areas are equal, and many apply some form of partial protection (Guidetti et al., 2014;Sala et al., 2012;Zupan et al., 2018). Such variation in regulations (e.g. gear and effort allowed) and enforcement is often associated with the response to protection (Edgar et al., 2014;Guidetti et al., 2014;Zupan et al., 2018). In practice, partial protection encompasses a wide variety of permitted use types (Zupan et al., 2018), as well as various levels of enforcement (Giakoumi et al., 2017;Guidetti et al., 2014;Sala et al., 2012). Here, we found high β-S values (i.e. spatial turnover due to rare and common species) in partially protected areas ( Figure S6). Although the data do not allow us to determine the underlying driver, these results suggest that siteto-site variation in permitted exploitation (i.e. regulation), enforcement or MPA effectiveness can increase the site-to-site variation in the fish community. This highlights the importance of quantifying site-to-site (i.e. spatial) variation in the response to protection, particularly where partial protection is a component of regional conservation efforts.
Regional conservation plans increasingly consider connectivity between individual protected areas as a key design feature (Green et al., 2015). The regional MPA system that we examined in this study represents an independently implemented ad hoc collection of MPAs, which were not established with a cohesive goal (see Grorud-Colvert et al., 2014). Adult movement for some of the exploited species in this study (i.e. Diplodus sargus, D. vulgaris, Epinephelus costae and E. marginatus) estimated using home ranges (i.e. the area where most time is spent foraging and resting), suggests that movement as adults is likely restricted to individual or, at most, adjacent sites within a given protected area (Di Franco et al., 2018). This means that protected sites in this study are likely relatively independent samples of adult populations. However, further empirical work will be required to examine whether the scale-dependent response of biodiversity to protection in networks designed with particular social or ecological goals (e.g. MPAs connected by larval, juvenile or adult dispersal) differs from those observed in this study.
Overall, our results show that analyses of multiple metrics across scales more fully reveals how biodiversity responds to protection. For the Mediterranean sites in this study, increased evenness played the predominant role in changes in biodiversity and site-to-site variation among the common species in the community was higher in protected areas. Identifying the drivers of these patterns will be an important next step for managing the Mediterranean MPAs in our study.
Additionally, the MPAs within the Mediterranean are typically small and cover a common pool of species. It would be revealing to examine whether our results hold across more heterogeneous MPA systems, that may consist of larger reserves, and cover multiple species pools and larger environmental gradients.