# Geographic distance, water circulation and environmental conditions shape the biodiversity of Mediterranean rocky coasts

Alex Rattray, Marco Andrello, Valentina Asnaghi, Stanislao Bevilacqua, Fabio Bulleri, Emma Cebrian, Mariachiara Chiantore, Joachim Claudet, Salud Deudero, Julian Evans, et al.

► To cite this version:

Alex Rattray, Marco Andrello, Valentina Asnaghi, Stanislao Bevilacqua, Fabio Bulleri, et al. "Geographic distance, water circulation and environmental conditions shape the biodiversity of Mediterranean rocky coasts." Marine Ecology Progress Series, 2016, 553, pp.1-11. 10.3354/meps11783. hal-01357935

HAL Id: hal-01357935

https://univ-perp.hal.science/hal-01357935

Submitted on 31 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
FEATURE ARTICLE

<figure><img src="image_2.png" /></figure>

# Geographic distance, water circulation and environmental conditions shape the biodiversity of Mediterranean rocky coasts

Alex Rattray<sup>1,2,*</sup>, Marco Andrello<sup>3</sup>, Valentina Asnaghi<sup>4</sup>, Stanislao Bevilacqua<sup>5</sup>, Fabio Bulleri<sup>1</sup>, Emma Cebrian<sup>6</sup>, Mariachiara Chiantore<sup>4</sup>, Joachim Claudet<sup>7,8</sup>, Salud Deudero<sup>9</sup>, Julian Evans<sup>10</sup>, Simonetta Fraschetti<sup>5</sup>, Giuseppe Guarnieri<sup>5</sup>, Luisa Mangialajo<sup>11,12</sup>, Patrick J. Schembri<sup>10</sup>, Antonio Terlizzi<sup>5,13</sup>, Lisandro Benedetti-Cecchi<sup>1</sup>

ABSTRACT: Ecological connectivity is important for effective marine planning and biodiversity conservation. Our aim was to identify factors important in influencing variation in benthic community structure on shallow rocky reefs in 2 regions of the Mediterranean Sea with contrasting oceanographic regimes. We assessed beta (β) diversity at 146 sites in the littoral and shallow sub-littoral from the Adriatic/ Ionian Seas (eastern region) and Ligurian/Tyrrhenian Seas (western region) using a null modelling approach to account for variation in species richness. The distance decay relationship between species turnover within each region and geographic distance by sea was determined using generalised linear models. Mantel tests were used to examine correlations between β diversity and connectivity by ocean currents, estimated from Lagrangian dispersal simulations. Variation in β diversity between sites was partitioned according to environmental and spatial components using a distance-based redundancy approach. Species turnover along a gradient of geographic distance was greater by a factor of 3 to 5 in the western region than the eastern region, suggesting lower connectivity between sites. β diversity was correlated with connectivity by ocean currents at both depths in the eastern region but not in the western region. The influ-

<figure><img src="image_3.png" /><figcaption>Diverse benthic assemblages on a shallow Mediterranean reef.

Photo: Lisandro Benedetti-Cecchi</figcaption></figure>

ence of spatial and environmental predictors of β diversity varied considerably between regions, but was similar between depths. Our results highlight the interaction of oceanographic, spatial and environmental processes influencing benthic marine β diversity. Persistent currents in the eastern region may be responsible for lower observed β diversity compared to the western region, where patterns of water circulation are more variable.

KEY WORDS: Beta diversity · Reefs · Spatial variation · MPA · Lagrangian modelling · Connectivity

<footnote n="1"><sup>1</sup>Department of Biology, University of Pisa, CoNISMa, Via Derna 1, 56126 Pisa, Italy</footnote>

<footnote n="2"><sup>2</sup>Centre of Integrative Ecology, School of Life and Environmental Sciences, Deakin University, Warrnambool, Victoria 3280, Australia</footnote>

<footnote n="3"><sup>3</sup>CEFE UMR 5175, CNRS - Université Paul-Valéry Montpellier - EPHE, Laboratoire Biogéographie et Écologie des Vertébrés, 1919 Route de Mende, 34293 Montpellier CEDEX 5, France</footnote>

<footnote n="4"><sup>4</sup>Department of Earth, Environment and Life Sciences, University of Genova, CoNISMa, Corso Europa 26, 16132 Genova, Italy</footnote>

<footnote n="5"><sup>5</sup>Department of Biological and Environmental Sciences and Technologies, University of Salento, CoNISMa, 73100 Lecce, Italy</footnote>

<footnote n="6"><sup>6</sup>Centre for Advanced Studies of Blanes (CEAB-CSIC), 17300 Blanes, Spain</footnote>

<footnote n="7"><sup>7</sup>National Center for Scientific Research, CRIOBE, USR 3278 CNRS-EPHE-UPVD, 66860 Perpignan, France</footnote>

<footnote n="8"><sup>8</sup>Laboratoire d'Excellence CORAIL, France</footnote>

<footnote n="9"><sup>9</sup>Instituto Español de Oceanografía, Centro Oceanográfico de Baleares, Moll de Ponent s/n, 07015 Palma de Mallorca, Spain</footnote>

<footnote n="10"><sup>10</sup>Department of Biology, University of Malta, Msida MSD 2080, Malta</footnote>

<footnote n="11"><sup>11</sup>Université Nice-Sophia Antipolis, EA 4228 ECOMERS, 28, avenue Valrose, 06108 Nice, France</footnote>

<footnote n="12"><sup>12</sup>Observatoire Océanologique de Villefranche sur mer, LOV CNRS UMR 7093, 06230 Villefranche-sur-mer, France</footnote>

<footnote n="13"><sup>13</sup>Stazione Zoologica Anton Dohrn, 80121 Napoli, Italy</footnote>

<footnote>*Corresponding author: alex.r@deakin.edu.au</footnote>
INTRODUCTION

Elucidating the mechanisms that control patterns
of diversity is a central objective of community ecol-
ogy (Hubbell 2001, Leibold et al. 2004), and has been
the focus of considerable research across terrestrial
(Tuomisto et al. 2003, De Cáceres et al. 2012), fresh-
water (Maloney & Munguia 2011, Angeler 2013) and
marine systems (Watson et al. 2012, Moritz et al.
2013). A better understanding of these mechanisms
is critically important for conservation planning,
which has historically focussed on prioritizing protec-
tion of areas containing high species richness and
threatened or endemic species, but often failed to
incorporate underlying variation in regional species
assembly driven by ecological connectivity, spatial
processes and environmental heterogeneity (Mc-
Knight et al. 2007).

Beta (β) diversity is a fundamental concept that addresses the spatial organisation of biodiversity. Although there are a multitude of approaches to defining and measuring β diversity, it is generally recognised as the change in species composition among sites in a geographic area of interest (Koleff et al. 2003, Legendre et al. 2005, Anderson et al. 2011). It can act as a measure of the change in species identities between 2 or more local assemblages, or the relationship between the diversity of the species pool for a given geographic region (γ diversity) and community diversity at local scales (α diversity) (Whittaker 1960). In each case, it reflects the range of processes acting on regional metacommunities that create and maintain patterns of biodiversity within them (de Juan et al. 2013). Anderson et al. (2011) identified 2 broad approaches to distinguishing β diversity: (1) *turnover* in species identities among sites along predefined spatial, environmental or temporal gradients or (2) *variation* in species identities among sites without reference to any specific gradient, both of which are useful for conservation planning (Legendre et al. 2005, Smith & Lundholm 2010, Legendre & De Cáceres 2013). Additionally, patterns of β diversity provide a measure of ecological connectivity that can be related to broader environmental, spatial and temporal gradients and used to untangle drivers behind patterns of community assembly (Borcard et al. 1992, Thrush et al. 2010).

In the Mediterranean region, as elsewhere in the
world, marine ecosystems are under mounting pres-
sure from the combined effects of environmental
stressors due to climate change and human impacts
(Bianchi & Morri 2000, Lejeusne et al. 2010). Efforts
to protect marine biodiversity through designation of

Marine Protected Areas (MPAs) have been under-
way for decades, although designation of MPAs has
often been based on opportunity rather than knowl-
edge of ecosystem functioning (Fraschetti et al. 2002),
leading to ad hoc systems rather than networks of
MPAs (Guilhaumon et al. 2015). The importance of
connected networks of reserves has been recognised
as a critical requirement to protecting the Mediter-
ranean marine ecosystem as a single interconnected
entity (Andrello et al. 2013, 2015, Lagabrielle et al.
2014), leading to internationally coordinated efforts
to plan for ecological connectivity between reserves
(Lubchenco et al. 2003, Grorud-Colvert et al. 2014).

In the present study, connectivity can be defined as
the exchange of individuals between populations
(Almany et al. 2009). It is widely recognised that high
connectivity increases resilience of individual MPAs
within a network by enhancing exogenous sources of
recruits, increasing genetic diversity and ultimately
maximising adaptive potential in the case of local
disturbance (Momigliano et al. 2015). Thus, con-
nected networks of MPAs can improve the ability of
ecosystems to face both natural and anthropogenic
perturbations. There is evidence that the potential of
connectivity generally decreases with geographic
distance (Cowen et al. 2006), but pathways and bar-
riers imposed by ocean currents, fronts and physico-
chemical conditions are also critical aspects to under-
standing the processes related to species' dispersal
(Shanks 2009). Most benthic species show benthic-
pelagic coupling of life cycles, where sedentary or
sessile adult stages rely on current-mediated advec-
tion for dispersal of pelagic juvenile forms (i.e. larvae
and propagules). Subsequent settlement success and
survival is dependent on the suitability of the receiv-
ing environment to species' requirements.

Assessing the relative contribution of distance,
ocean currents and suitability of environmental con-
ditions to patterns of β diversity represents an impor-
tant prerequisite for MPA network design (Watson
et al. 2011b, Rassweiler et al. 2012). Planning con-
nected networks of MPAs requires an understanding
of how each of these factors contributes to network
connectivity as a whole. While a number of studies
have explored the interrelated influence of these fac-
tors at the population level, their role in community-
level dynamics is less well understood. This is partly
due to the difficulties of incorporating the wide vari-
ation in species life histories and dispersal traits into
a single conceptual model (Moritz et al. 2013).

Recent refinement of ocean circulation models and
improvements in computing power have allowed for a
better understanding of the degree of hydrodynamic
connectivity between populations (Andrello et al. 2013, 2015) and communities (Watson et al. 2011a). Similarly, ocean color satellite imagery provides synoptic information on the distribution of biologically relevant environmental conditions across temporal and geographic scales not readily attainable by direct sampling alone. In this study, we combined this information to examine the influence of geographic distance, current connectivity potential and environmental variability on patterns of β diversity in 2 regions of the Mediterranean Sea with contrasting oceanographic regimes. Persistent oceanographic features in the southern Adriatic and Ionian Seas and more variable, seasonal patterns in the northern Tyrrhenian and Ligurian Seas provide a unique opportunity to assess the effects of connectivity potential on β diversity. We specifically addressed the following questions: (1) Do patterns of β diversity vary consistently between depths and regions of the Mediterranean Sea? (2) Which measure of connectivity, geographic or connection by ocean currents, best correlates with β diversity? (3) What environmental factors determine the variation in β diversity between depths and regions, and what are their relative contributions?

# MATERIALS AND METHODS

## Biological data

We analyzed β diversity in assemblages of algae and sessile invertebrates from rocky substrates at 2 depths on coastlines in the Ligurian, Tyrrhenian, Adriatic, Ionian and Central Mediterranean Seas (Fig. 1). Data for littoral areas (0 to -0.2 m relative to mean low water level) were obtained for 60 replicated sites (5100 quadrats) from surveys carried out between 2000 and 2008. A total of 86 sites (4150 quadrats) in the shallow sublittoral (5 to 7 m) were sampled between 1999 and 2011. At littoral sites, visual estimates of species cover were made directly in the field using 20 × 20 cm (0.04 m²) quadrats. For the sublittoral, an area of 16 × 25 cm (0.04 m²) was photographically sampled, and then species cover was visually estimated from slides under magnification. For comparative regional analyses, sites were assigned to 2 regions. The western region (35 littoral and 48 sublittoral sites) encompassed the eastern Ligurian and northern Tyrrhenian Seas including MPAs at Portofino, Cinque Terre, and the Tuscan

<figure><img src="image_4.png" /><figcaption>Fig. 1. Littoral and sublittoral study locations for analysis of biodiversity distribution. Western and eastern sites used in regional comparisons are defined by red bounding boxes. Sites in the northern Adriatic and central Mediterranean (Malta) were included in a global analysis of species turnover</figcaption></figure>
Archipelago (Capraia, Montecristo, Giannutri and Elba Islands). The eastern region (25 littoral and 38 sublittoral sites) included the southern Adriatic and northern Ionian Seas including MPAs at Porto Cesareo, Torre Guaceto, and the Tremiti Islands. An additional 10 sites in Malta and the northern Adriatic (Miramare MPA) were not assigned to either of the regions but were included in a global analysis of β diversity combining sites from both regions.

## β diversity as variation

Observed β diversity (β₀) was derived separately for all pairs of littoral and sublittoral sites in each region using Jaccard's coefficient of dissimilarity (dⱼ):

$$d_J = \frac{(b+c)}{(a+b+c)} \quad (1)$$

in which a is the number of species shared by 2 sites, b is the number of species in the other site, but not in the focal one, and c is the number of species unique to the focal site. Estimates of β diversity are concomitantly dependent on variation in local species richness (α diversity) (Chase et al. 2011, Kraft et al. 2011). To separate underlying assembly processes that create β diversity from the influence of α diversity, a null modelling approach was used. A 'null community' for each site was derived by randomly sampling from the full species pool within each depth category and region. Sample size was equal to the number of taxa observed at a given site, and the probability of each taxon being sampled was weighted proportionally to its relative cover. Jaccard dissimilarities were calculated and averaged between all possible pairs of sites over 9999 iterations of the resampling procedure to obtain a null estimate of between-site β diversity (βₙ). Resulting differences (βₚₐᵥ) between observed and expected beta diversity under the null expectation (β₀ – βₙ) were used in all subsequent analyses.

## Environmental covariates

A suite of satellite-derived oceanographic variables were selected as environmental predictors of β diversity. Level 3 Moderate Resolution Imaging Spectro-radiometer (MODIS) data with spatial resolution of ~4 km at the equator and with daily temporal resolution were acquired from NASA's Goddard Space Flight Centre (http://oceancolor.gsfc.nasa.gov/). Cumulative climatologies between the years 2000 (the earliest available MODIS data) and 2011 were

derived for night-time sea surface temperature (SST), sea surface salinity (Salinity), net primary productivity (NPP) (Behrenfeld & Falkowski 1997), diffuse light attenuation coefficient at 490 nm (kₐ₄₉₀), Chlorophyll a concentration (CHL) and concentration of colored dissolved organic matter (CDOM). For each of the satellite products the mean, standard deviation around the mean, mean annual range and highest monthly mean were calculated over the full time-series for each pixel corresponding to a single sampling site. A long-term averaging approach was taken in order to accommodate the differences in field sampling times. Sites beyond the extent of the satellite climatologies (n = 7) were assigned the value of the nearest adjacent pixel without interpolation. Prior to analysis, multi-collinearity and correlation of environmental predictor variables was assessed for each region using a step-wise procedure, where at each iteration the variable with the greatest variance inflation factor (VIF) was removed until remaining covariates (SSTₐugust, SSTₐnn.Ranger, CDOMₚMean, CDOMₚSD, kₐ₄₉₀ₚSD) displayed VIF values <5. Since sampling sites were located within and outside MPAs, we included protection status (fully protected vs. partially protected and unprotected) as an environmental covariate in the analyses to account for any effects of protection on β diversity.

## Geographic distance

Pair-wise distances were calculated between sites for each region. Distances by sea between each site were calculated using a least-cost paths approach (hereafter least-cost distance). Land areas were masked using the European Environment Agency coastline polygon (1:100 000) (www.eea.europa.eu/legal/copyright) and defined as a barrier in distance calculations. Given the complex geometry of coastlines in the study regions this was considered a more realistic estimate of effective spatial separation than the more commonly used Great Circle or Euclidean distances.

## Connectivity by ocean currents

A current connectivity matrix was derived using Lagrangian particle simulations in order to estimate potential connectivity by ocean currents. Three-dimensional sea current velocities and directions were obtained from the 1/12th degree (6 to 8 km resolution) NEMOMED12 model (https://hal-ensta.archives-
ouvertes.fr/hal-01138958/). Due to the inherent difficulties in developing a model representing the multitude of species' dispersal strategies, no weighting or biotic parameterization was applied to the model. Simulation of Lagrangian particle transit times therefore represents the movement of particles of water between sites. In this vein, current connectivity was defined as the minimum particle transit time between pairs of sites. A detailed description of the hydrodynamic simulations can be found in Andrello et al. (2013).

## Spatial analysis

Spatial predictors were generated as Moran's eigenvector maps (MEMs) (Dray et al. 2006) by a principal coordinate analysis (PCoA) on the truncated spatial distance matrix. Eigenvectors represent spatial relationships among sampling sites across a range of spatial scales. Primary eigenvectors define relationships between sampling sites over broad spatial scales, and subsequent eigenvectors provide information at progressively finer scales. Additionally, the decomposition of the spatial relationships between sampling locations was weighted using the current connectivity matrix (days) described above. The spatial predictors are therefore representative of site connectivity by currents across multiple scales.

## β diversity as turnover

Distance decay of similarity is a measure describing the contribution of geographic distance to changes in community similarity. Based on the distance decay relationship, measures of initial similarity ($S_0$) and halving distance ($D_h$) can be calculated. $S_0$ reflects β diversity at small spatial extents, whereby higher values of $S_0$ indicate lower β diversity at small scales. $D_h$ indicates the distance at which $S_0$ halves, and provides relevant information regarding the spatial scale of species turnover. Pair-wise Jaccard similarities ($S_{ij} = 1 - d_{ij}$) (Anderson et al. 2013) were plotted against least-cost distance with distance decay relationships estimated by generalised linear models (GLMs) with binomial observation error and a log-link function (Millar et al. 2011), taking the form:

$$ \log(E[S]) = \alpha - \beta D \quad (2) $$

where $\log(E[S])$ denotes the expected similarity, $D$ is least-cost distance (km), and $\alpha$ and model slope, $\beta$, are the parameters of the log-linear model. As pairwise similarities are not independent, parameter esti-

mates and their standard errors for $\alpha$, $\beta$ and $D_h$ were obtained by a leave-one-out jackknife procedure using 9999 permutations (Millar et al. 2011). Because species turnover was positively related to protection status (Table 2), we performed a sensitivity analysis to evaluate the potential confounding effect due to having a larger number of sites within MPAs in the western than in the eastern region. Parameters $\alpha$, $\beta$ and $D_h$ were estimated for residuals of a regression of Jaccard similarity against level of protection for depths within regions. Since the results and their interpretation were not substantially different from the results obtained using the approach outlined above, they have not been included in the main text (See Table S1 in the Supplement at www.int-res.com/articles/suppl/m553p001_supp.pdf). Distance decay relationships were similarly derived between pairwise Jaccard similarity and travel time (minimum and mean) estimated from Lagrangian particle simulations and are presented in Fig. S1 in the Supplement.

## Correlation of species turnover with current connectivity and least-cost distance

To assess the relative influence of connectivity by ocean currents and geographic distance on variation in β diversity, Mantel tests were used to test the significance of the Pearson product-moment coefficient (r) between pair-wise Jaccard dissimilarity ($d_J$) and matrices of environmental distance, least-cost distance and current connectivity with significance assessed using 9999 random permutations. The Euclidean environmental distance matrix was calculated between sites using non-collinear (VIF < 5) environmental covariates ($SST_{August}$, $SST_{Ann.Range}$, $CDOM_{Mean}$, $CDOM_{SD}$, $K_d490_{SD}$). Complementary partial Mantel tests were then done to determine correlations between β diversity, least-cost and current distance whilst controlling for the effect of environmental distance.

## Partitioning the contribution of environmental and spatial predictors

A variation partitioning approach (Borcard et al. 1992, Legendre et al. 2005) was used to determine the relative contribution of each set of environmental and spatial predictors (i.e. MEMs) to β diversity. This was carried out by means of 3 distance-based redundancy analyses (dbRDA) considering both spatial and environmental predictors together or separately.
The total percentage of variation explained by the model (adjusted R²) was partitioned into unique and common contributions of each of the sets of environmental and spatial predictors. Fractions of variation describing the independent effects of spatial and environmental processes (along with their non-separable fractions and unexplained variation) were then obtained by subtraction (e.g. Borcard et al. 1992, Peres-Neto et al. 2006).

## Environmental drivers of β diversity

Finally, we used multiple regression to determine the contribution of individual environmental covariates to predicting β diversity for combined littoral sites from both regions (including those in the northern Adriatic and Malta), and combined sublittoral sites (including those in the northern Adriatic). In order to reduce pairwise dissimilarities to a single measure of β diversity for each site, the average Jaccard dissimilarity was calculated among all pairs of sites (Anderson et al. 2013, de Juan et al. 2013). That is, the β diversity for a given focal site was computed as the mean of all the pair-wise dissimilarities between the focal site and each of the other sites. As

significant positive spatial autocorrelation (α < 0.05) was present in the residuals of the ordinary least squares (OLS) regression, a spatial generalised least squares (GLS) regression was used. The GLS is derived from standard regression, but differs in the incorporation of the spatial structure directly into the residuals, assuming in this study an exponential relationship between the error term and geographic distance (Diniz Filho & Bini 2005).

All statistical analyses were done using R version 3.1.2 (R Development Core Team 2014) (packages 'vegan', 'nlme', 'spacemakeR', 'simba', 'spam') with code obtained from Millar et al. (2011) and Chase et al. (2011).

# RESULTS

## β diversity as turnover

Initial similarities were relatively high between littoral sites in both the eastern (0.239) and western (0.293) regions, and low between sublittoral sites (0.161 and 0.120 in the eastern and western regions respectively). This pattern of lower β diversity at small spatial extents in littoral habitats was reflected

<figure><img src="image_5.png" /><figcaption>Fig. 2. Decay in site similarities (Jaccard) of benthic communities with increasing distance in (A) littoral and (B) sublittoral sampling sites in the western region; and (C) littoral and (D) sublittoral sampling sites in the eastern region. The fitted models from the binomial generalised linear models (GLMs) with log-link are shown with corresponding parameter estimates (β) and halving distances ($D_h$).</figcaption></figure>
in broader-scale species turnover, although there was more variation between regions (Fig. 2). The rate of species turnover (i.e. slope β) with distance was greater for the western region by a factor of ~3 between littoral sites and by a factor of ~5 for sublittoral sites. Corresponding D<sub>h</sub> (the mean distance by sea at which community similarity halves) varied between depths and regions. These were estimated at 475 ± 168.9 and 1344.5 ± 592.8 km for littoral sites, and 138.8 ± 10.5 and 639.2 ± 157.1 km for sublittoral sites in the western and eastern study regions respectively.

## Correlation of species turnover with current connectivity and least-cost distance

Pair-wise β diversity between littoral sites in the western region was strongly positively correlated with least-cost distance (r = 0.670, p < 0.001) and shared weaker, but significant correlations with current distance (r = 0.185, p < 0.001) and environmental distance (r = 0.226, p < 0.01) (Table 1). Least-cost distance was likewise strongly correlated with β diversity between sublittoral sites in this region (r = 0.675, p < 0.001), as was environmental distance (r = 0.675, p < 0.001). In the eastern region, β diversity was weakly but significantly correlated with current distance (r = 0.248, p < 0.001), but not with either distance by sea or environmental distance.

<table><caption>Table 1. Correlations between β diversity (i.e. community variation), ocean distance, environmental distance and current connectivity between pairs of sampling sites and partial correlations after controlling for the effects of environmental distance. The statistical significance of comparisons was assessed using Mantel and partial Mantel tests based on Pearson's product-moment correlation using 9999 permutations; \*p < 0.05, \*\*p < 0.01</caption><thead><tr><th></th><th>N pairs</th><th>Least-cost<br>distance (km)</th><th>Current distance<br>(days)</th><th>Environmental<br>distance</th></tr></thead><tbody><tr><td colspan="5"><strong>Correlations</strong></td></tr><tr><td colspan="5"><strong>Littoral</strong></td></tr><tr><td>West</td><td>612</td><td>0.670**</td><td>0.185*</td><td>0.226**</td></tr><tr><td>East</td><td>312</td><td>0.161</td><td>0.248**</td><td>0.122</td></tr><tr><td colspan="5"><strong>Sublittoral</strong></td></tr><tr><td>West</td><td>968</td><td>0.675**</td><td>0.097</td><td>0.567*</td></tr><tr><td>East</td><td>882</td><td>0.135</td><td>0.142**</td><td>0.1232*</td></tr><tr><td colspan="5"><strong>Partial correlations, corrected for environmental distance</strong></td></tr><tr><td colspan="5"><strong>Littoral</strong></td></tr><tr><td>West</td><td>612</td><td>0.661**</td><td>0.168*</td><td>-</td></tr><tr><td>East</td><td>312</td><td>0.109</td><td>0.217**</td><td>-</td></tr><tr><td colspan="5"><strong>Sublittoral</strong></td></tr><tr><td>West</td><td>968</td><td>0.455**</td><td>-0.048</td><td>-</td></tr><tr><td>East</td><td>882</td><td>0.083</td><td>0.112*</td><td>-</td></tr></tbody></table>

## Partitioning the contribution of environmental and spatial predictors

The relative influence of spatial and environmental processes in describing community structure varied substantially between study regions, but was remarkably similar between littoral and sublittoral sites within each region (Fig. 3). In the western region, spatial and environmental predictors accounted for approximately equal portions of variation in community structure (7 to 10%) with large overlap between the two (51 to 53%). In contrast, spatial gradients mediated by current connectivity were significantly more important than environmental predictors in describing community variation in the eastern region for both littoral and sublittoral sites (19% in both cases). Environmental covariates alone explained only 5% of variation in β diversity between littoral sites and none of the variation between sublittoral sites.

## Environmental drivers of β diversity

Multiple regression was conducted to test whether satellite-derived environmental covariates predicted species turnover for all littoral (including those in the northern Adriatic and Malta), and all sublittoral sites (including those in the northern Adriatic). Based on Akaike's information criterion

(AIC), multivariate models based on spatial GLS regression showed a better description of environmental variables on species turnover than equivalent OLS models, which do not take spatial autocorrelation into account (Table 2). For littoral sites, SST variables (SST<sub>Ann.Range</sub>, SST<sub>Mean</sub>) and protection status were significant predictors of species turnover (F<sub>6,60</sub> = 7.19, p < 0.001, R<sup>2</sup><sub>adj</sub> = 0.36). In the sublittoral sites, long-term variation of CDOM (CDOM<sub>SD</sub>) and mean annual range of SST (SST<sub>Ann.Range</sub>) were the best environmental correlates of species turnover (F<sub>6,78</sub> = 7.45, p < 0.0001, R<sup>2</sup><sub>adj</sub> = 0.33).
<figure><img src="image_6.png" /><figcaption>Fig. 3. Variation partitioning (Venn diagrams) of unique contribution of spatial (Sp.) and environmental (Env.) components to benthic β diversity distribution, for the littoral and sublittoral in the western and eastern regions. Overlapping fractions represent the shared variation between the environmental and spatial components and the residual variation unexplained by the multivariate model. p-values showing the significance of each fraction of variation were estimated with 199 permutations under the full model (*p < 0.05)</figcaption></figure>

DISCUSSION

The present study provides an example of how decoupling the influence of ocean currents, spatial and environmental drivers can provide important insights into the factors affecting regional patterns of β diversity. We found that regional scale (100s of km) patterns in species turnover of rocky reef benthos in the Mediterranean differ between regions and depths, the relationship between species turnover and connectivity by ocean currents was regionally variable, and finally, the relative influence of spatial and environmental processes in describing community structure varied substantially between study regions, but was remarkably similar between littoral and sublittoral sites within each region.

Mechanisms controlling the decay of community similarity with distance have been broadly, but not exclusively, attributed to 3 key factors (Nekola & White 1999, Soininen et al. 2007): niche-based processes, spatial configuration of the environment and stochastic processes posited by neutral theory (e.g. dispersal limitation) (Hubbell 2001). In agreement with previous studies (Watson et al. 2011a, Moritz et al. 2013), our findings provide support for each of these factors, varying in degree by region and across depths. We found that decay of community similarity with distance was comparatively lower for the southern Adriatic and Ionian Seas, suggesting greater regional community homogeneity than in the Tyr-

<table><caption>Table 2. Partial regression coefficients of the multiple regression models (b), standard errors of the estimates (SE), t-statistics and associated p-values for Mediterranean benthic β diversity at 2 depths regressed against satellite-derived environmental variables. An exponential model was assumed in generalised least squares (GLS) regression to describe spatial patterns in model residuals. Akaike's information criteria (AIC) values are reported from spatial GLS models (AIC<sub>GLS</sub>) and non-spatial ordinary least squares (OLS) models (AIC<sub>OLS</sub>) containing the same variables; adjusted R<sup>2</sup> (R<sup>2</sup><sub>adj</sub>) are reported from non-spatial OLS models. For the littoral zone: AIC<sub>GLS</sub> = -212.2, AIC<sub>GLS</sub> = -244.6; R<sup>2</sup> = 0.42, R<sup>2</sup><sub>adj</sub> = 0.36. For the sublittoral zone: AIC<sub>OLS</sub> = -322.1, AIC<sub>GLS</sub> = -352.7; R<sup>2</sup> = 0.36, R<sup>2</sup><sub>adj</sub> = 0.33. SST: sea surface temperature; CDOM: colored dissolved organic matter; k<sub>d</sub>490: diffuse attenuation coefficient for downwelling irradiance at 490 nm; 'other' protection status: partially protected and unprotected. *p < 0.05; **p < 0.01</caption><thead><tr><th rowspan="2">Variable</th><th colspan="4">Littoral (0–0.2 m)</th><th colspan="4">Sublittoral (5–7 m)</th></tr><tr><th>b</th><th>SE</th><th>t</th><th>p</th><th>b</th><th>SE</th><th>t</th><th>p</th></tr></thead><tbody><tr><td>Intercept</td><td>-1.06</td><td>0.373</td><td>-2.86</td><td>0.006**</td><td>0.159</td><td>0.274</td><td>0.58</td><td>0.563</td></tr><tr><td>SST</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean</td><td>0.047</td><td>0.016</td><td>2.85</td><td>0.006**</td><td>-0.009</td><td>0.015</td><td>-0.59</td><td>0.556</td></tr><tr><td>Annual range</td><td>0.024</td><td>0.009</td><td>2.74</td><td>0.008**</td><td>0.011</td><td>0.005</td><td>2.43</td><td>0.017*</td></tr><tr><td>CDOM</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean</td><td>0.013</td><td>0.008</td><td>1.59</td><td>0.117</td><td>0.001</td><td>0.007</td><td>0.16</td><td>0.873</td></tr><tr><td>SD</td><td>0.006</td><td>0.014</td><td>-0.42</td><td>0.675</td><td>-0.029</td><td>0.012</td><td>-2.42</td><td>0.019*</td></tr><tr><td>k<sub>d</sub>490</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SD</td><td>0.297</td><td>0.525</td><td>0.56</td><td>0.573</td><td>0.621</td><td>0.381</td><td>1.63</td><td>0.108</td></tr><tr><td>Protection status</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fully protected vs. other</td><td>0.028</td><td>0.009</td><td>2.84</td><td>0.006**</td><td>-0.009</td><td>0.007</td><td>-1.19</td><td>0.239</td></tr></tbody></table>
rhenian/Ligurian region. Halving distances in the eastern region were also large relative to the spatial scale of the study, suggesting a homogeneous and potentially well-connected system. The shallower rates of decay observed in the eastern region may be attributable to lower environmental variability, conforming to the concept of species sorting by environmental niche processes. This contention was further supported by regional differences in the contribution of environmental variables to partitioned variation in community similarity.

Spatial processes explained the greatest fraction of variation of β diversity in the eastern region. While it is difficult to relate spatial structures to the ecological processes that create them (Hubbell 2001, Smith & Lundholm 2010), significant correlations were also found between pair-wise measures of β diversity along a gradient of connectivity by currents. This may indicate that dispersal potential is an important factor determining community organisation in the region. The opposite is true for variation among sites in the western region. Here, we observed a strong and significant effect of geographic distance on decreasing similarity between communities, but only a relatively minor effect of current connectivity in littoral communities and no significant effect in sublittoral communities. While this may result from unmeasured anthropic disturbances such as localised extractive activities (Fraschetti et al. 2001) and urbanisation effects (Benedetti-Cecchi et al. 2001), at least some of these discrepancies between regions may be explained by differences in oceanographic regimes. The southern Adriatic (eastern region) is characterised by a generally stable cyclonic cell, largely controlled by the topography of the South Adriatic Pit, between the western and eastern Adriatic boundary currents (Poulain 2001). In contrast, current patterns along the Tyrrhenian and Ligurian coasts (western region) display a generally weak and variable northward trend as the Eastern Corsican Current moves through the Corsican channel into the Liguro-Provençal Basin (Astraldi et al. 1990). Several authors have noted an anti-cyclonic eddy persistently forming in the summer months at the latitude of the island of Capraia, which may seasonally limit or even reverse transport through the Corsican Channel (Astraldi et al. 1990, Iacono et al. 2013). Andrello et al. (2013) reported only north-south directional connectivity between the 2 basins in Lagrangian simulations between MPAs. The weaker and more variable current regime connecting the northern Tyrrhenian and southern Ligurian seas, coupled with complex seasonal wind-driven dynamics, may therefore consti-

tute a barrier to dispersal between northern and southern sites in the region. An important caveat to these findings is that our biophysical model was not calibrated according to the specific life histories of any single species or group of species. Species found within MPAs display a multitude of dispersal strategies which are dependent on many factors, including propagule morphology, movement characteristics, pelagic duration and timing of release (Shanks 2009), that may relate to habitat features on a different scale (Thrush et al. 2011). A number of theories have been proposed to explain dynamics and scales of metacommunity structure (Leibold et al. 2004). Deterministic and stochastic processes responsible for metacommunity assembly may interact over a range of spatial and temporal scales. In generally homogeneous environments, dispersal may be the more important organisational mechanism, while characteristics of the receiving environment may predominate where conditions show small-scale variability. The available evidence suggests that this is the case in the present study. Despite their similar latitudinal extents, environmental variation was lower in the eastern region than in the western region and was largely unrelated to patterns of species turnover, but significantly correlated with patterns in ocean currents.

Environmental variation as estimated from long-term satellite imagery indicated that sublittoral species turnover was best predicted by a negative relationship with CDOM variability and a positive relationship with annual variation in SST. Similarly, SST (mean and annual range) shared a significant positive relationship with species turnover at littoral sites, when considered against all other variables. Temperature is a well-documented correlate of broad-scale marine diversity gradients, although its effects may vary over finer scales according to local variation in primary productivity and community trophodynamics (Roy et al. 1998). Concentrations of CDOM in coastal waters vary according to *in situ* creation of fulvic acids (Hulatt et al. 2009), or to presence of humic acids in terrigenous effluent. Higher concentrations of CDOM dominate the absorption of both visible and blue spectrum UV light in the water column, potentially limiting photosynthetically available radiation to the sublittoral benthos. Modification of the underwater light regime can impact biomass and species richness of macroalgae (Desmond et al. 2015), which act as foundation species, providing habitat and hosting complex ecosystem processes to a suite of associated species.

A better understanding of the processes underlying patterns of regional biodiversity is highly relevant to
marine spatial planning, especially when assessing
candidate sites for protection and inclusion in net-
works of MPAs. A significant finding of this study is
that differences in patterns of community similarity
do not appear to be driven uniformly by the same
processes. Instead, the measured drivers behind vari-
ation in community structure were found to differ be-
tween regions, and to a lesser extent between depths.
Current connectivity, geographic distance and envi-
ronmental variability were all linked to patterns in β
diversity at varying degrees of importance across re-
gions and depths. This implies that strategies ac-
counting for this wide range of connectivity drivers
need to be considered when planning a network of
MPAs aimed at ensuring connectivity between units.
In the example of benthic communities in the study
systems considered here, higher priority would need
to be given to more heterogeneous regions, where
species turnover is greater over smaller spatial scales.
That is, units in the western region would need to
be closer together to capture a similar gradient of
change in diversity relative to those in the eastern re-
gion. Additionally, connectivity potential of candidate
habitats for protection is an important consideration
in network spatial arrangement. Our results suggest
that sublittoral community structure changes more
rapidly over distance than it does in the littoral.
Therefore, evaluation of unit connectivity should be
predicated on habitats that display the greatest com-
munity heterogeneity through space.

Acknowledgements. This work was supported by the European Community's 7th Framework Programme (FP7/2007–2013) under Grant Agreement No. 287844 for the project 'Towards COast to COast NETworks of marine protected areas (from the shore to the high and deep sea), coupled with sea-based wind energy potential (CoCoNet)'. The funding body had no role in study design, data collection and analysis, or preparation of the manuscript. Author contributions: A.R. and L.B.C. conceived the work and conducted the analysis; M.A. conducted the biophysical modeling; A.R. led the writing of the manuscript. All authors contributed data and reviewed the manuscript. We acknowledge and appreciate the effort of our reviewers for their invaluable assistance in improving the manuscript.

LITERATURE CITED

* Almany GR, Connolly SR, Heath DD, Hogan JD and others (2009) Connectivity, biodiversity conservation and the design of marine reserve networks for coral reefs. Coral Reefs 28:339–351

* Anderson MJ, Crist TO, Chase JM, Vellend M and others (2011) Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist. Ecol Lett 14:19–28

* Anderson MJ, Tolimieri N, Millar RB (2013) Beta diversity of

demersal fish assemblages in the north-eastern Pacific: interactions of latitude and depth. PLoS ONE 8:e57918

* Andrello M, Mouillot D, Beuvier J, Albouy C, Thuiller W, Manel S (2013) Low connectivity between Mediterranean Marine Protected Areas: a biophysical modeling approach for the dusky grouper *Epinephelus marginatus*. PLoS ONE 8:e68564

* Andrello M, Jacobi MN, Manel S, Thuiller W, Mouillot D (2015) Extending networks of protected areas to optimize connectivity and population growth rate. Ecography 38: 273–282

* Angeler DG (2013) Revealing a conservation challenge through partitioned long-term beta diversity: increasing turnover and decreasing nestedness of boreal lake meta-communities. Divers Distrib 19:772–781

* Astraldi M, Gasparini GP, Manzella GMR, Hopkins TS (1990) Temporal variability of currents in the eastern Ligurian Sea. J Geophys Res 95:1515–1522

* Behrenfeld MJ, Falkowski PG (1997) Photosynthetic rates derived from satellite-based chlorophyll concentration. Limnol Oceanogr 42:1–20

* Benedetti-Cecchi L, Pannacciulli F, Bulleri F, Moschella P, Airoldi L, Relini G, Cinelli F (2001) Predicting the consequences of anthropogenic disturbance: large-scale effects of loss of canopy algae on rocky shores. Mar Ecol Prog Ser 214:137–150

* Bianchi CN, Morri C (2000) Marine biodiversity of the Mediterranean Sea: situation, problems and prospects for future research. Mar Pollut Bull 40:367–376

* Borcard D, Legendre P, Drapeau P (1992) Partialling out the spatial component of ecological variation. Ecology 73: 1045–1055

* Chase JM, Kraft NJB, Smith KG, Vellend M, Inouye BD (2011) Using null models to disentangle variation in community dissimilarity from variation in α-diversity. Ecosphere 2:art24

* Cowen RK, Paris CB, Srinivasan A (2006) Scaling of connectivity in marine populations. Science 311:522–527

* De Cáceres M, Legendre P, Valencia R, Cao M and others (2012) The variation of tree beta diversity across a global network of forest plots. Glob Ecol Biogeogr 21: 1191–1202

* de Juan S, Thrush SF, Hewitt JE (2013) Counting on beta-diversity to safeguard the resilience of estuaries. PLoS ONE 8:e65575

* Desmond MJ, Pritchard DW, Hepburn CD (2015) Light limitation within southern New Zealand kelp forest communities. PLoS ONE 10:e0123676

Diniz Filho JAF, Bini LM (2005) Modelling geographical patterns in species richness using eigenvector based spatial filters. Glob Ecol Biogeogr 14:177–185

* Dray S, Legendre P, Peres-Neto PR (2006) Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol Modell 196: 483–493

* Fraschetti S, Bianchi CN, Terlizzi A, Fanelli G, Morri C, Boero F (2001) Spatial variability and human disturbance in shallow subtidal hard substrate assemblages: a regional approach. Mar Ecol Prog Ser 212:1–12

* Fraschetti S, Terlizzi A, Micheli F, Benedetti-Cecchi L, Boero F (2002) Marine Protected Areas in the Mediterranean Sea: objectives, effectiveness and monitoring. Mar Ecol (Berl) 23:190–200

* Grorud-Colvert K, Claudet J, Tissot BN, Caselle JE and others (2014) Marine Protected Area networks: assessing
whether the whole is greater than the sum of its parts. PLoS ONE 9:e102298

► Guilhaumon F, Albouy C, Claudet J, Velez L and others (2015) Representing taxonomic, phylogenetic and func- tional diversity: new challenges for Mediterranean marine-protected areas. Divers Distrib 21:175–187

Hubbell SP (2001) The unified neutral theory of biodiversity and biogeography. Princeton University Press, Princeton, NJ

► Hulatt CJ, Thomas DN, Bowers DG, Norman L, Zhang C (2009) Exudation and decomposition of chromophoric dissolved organic matter (CDOM) from some temperate macroalgae. Estuar Coast Shelf Sci 84:147–153

► Iacono R, Napolitano E, Marullo S, Artale V, Vetrano A (2013) Seasonal variability of the Tyrrhenian Sea surface geostrophic circulation as assessed by altimeter data. J Phys Oceanogr 43:1710–1732

► Koleff P, Gaston KJ, Lennon JJ (2003) Measuring beta diver- sity for presence–absence data. J Anim Ecol 72:367–382

► Kraft NJB, Comita LS, Chase JM, Sanders NJ and others (2011) Disentangling the drivers of β diversity along latitudinal and elevational gradients. Science 333: 1755–1758

► Lagabrielle E, Crochelet E, Andrello M, Schill SR, Arnaud- Haond S, Alloncle N, Ponge B (2014) Connecting MPAs – eight challenges for science and management. Aquat Conserv Mar Freshw Ecosyst 24:94–110

► Legendre P, De Cáceres M (2013) Beta diversity as the vari- ance of community data: dissimilarity coefficients and partitioning. Ecol Lett 16:951–963

► Legendre P, Borcard D, Peres-Neto PR (2005) Analyzing beta diversity: partitioning the spatial variation of com- munity composition data. Ecol Monogr 75:435–450

► Leibold MA, Holyoak M, Mouquet N, Amarasekare P and others (2004) The metacommunity concept: a framework for multi-scale community ecology. Ecol Lett 7:601–613

Lejeusne C, Chevaldonne P, Pergent-Martini C, Boudou- resque CF, Perez T (2010) Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea. Trends Ecol Evol 25:250–260

► Lubchenco J, Palumbi SR, Gaines SD, Andelman S (2003) Plugging a hole in the ocean: the emerging science of marine reserves. Ecol Appl 13:3–7

► Maloney KO, Munguia P (2011) Distance decay of similarity in temperate aquatic communities: effects of environ- mental transition zones, distance measure, and life histo- ries. Ecography 34:287–295

► McKnight MW, White PS, McDonald RI, Lamoreux JF, Sechrest W, Ridgely RS, Stuart SN (2007) Putting beta- diversity on the map: broad-scale congruence and coin- cidence in the extremes. PLoS Biol 5:e272

► Millar RB, Anderson MJ, Tolimieri N (2011) Much ado about nothings: using zero similarity points in distance-decay curves. Ecology 92:1717–1722

► Momigliano P, Harcourt R, Stow A (2015) Conserving coral reef organisms that lack larval dispersal: Are networks of

Marine Protected Areas good enough? Front Mar Sci 2: 1–5

Moritz C, Meynard CN, Devictor V, Guizien K, Labrune C, Guarini JM, Mouquet N (2013) Disentangling the role of connectivity, environmental filtering, and spatial struc- ture on metacommunity dynamics. Oikos 122:1401–1410

► Nekola JC, White PS (1999) The distance decay of similarity in biogeography and ecology. J Biogeogr 26:867–878

► Peres-Neto PR, Legendre P, Dray S, Borcard D (2006) Varia- tion partitioning of species data matrices: estimation and comparison of fractions. Ecology 87:2614–2625

► Poulain PM (2001) Adriatic Sea surface circulation as derived from drifter data between 1990 and 1999. J Mar Syst 29:3–32

R Development Core Team (2014) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna

► Rassweiler A, Costello C, Siegel DA (2012) Marine protected areas and the value of spatially optimized fishery man- agement. Proc Natl Acad Sci USA 109:11884–11889

► Roy K, Jablonski D, Valentine JW, Rosenberg G (1998) Marine latitudinal diversity gradients: tests of causal hypotheses. Proc Natl Acad Sci USA 95:3699–3702

► Shanks AL (2009) Pelagic larval duration and dispersal distance revisited. Biol Bull (Woods Hole) 216:373–385

► Smith TW, Lundholm JT (2010) Variation partitioning as a tool to distinguish between niche and neutral processes. Ecography 33:648–655

► Soininen J, McDonald R, Hillebrand H (2007) The distance decay of similarity in ecological communities. Ecography 30:3–12

► Thrush SF, Hewitt JE, Cummings VJ, Norkko A, Chiantore M (2010) β-diversity and species accumulation in Antarc- tic coastal benthos: influence of habitat, distance and productivity on ecological connectivity. PLoS ONE 5: e11899

► Thrush SF, Chiantore M, Asnaghi V, Hewitt J, Fiorentino D, Cattaneo-Vietti R (2011) Habitat–diversity relationships in rocky shore algal turf infaunal communities. Mar Ecol Prog Ser 424:119–132

► Tuomisto H, Ruokolainen K, Yli-Halla M (2003) Dispersal, environment, and floristic variation of western Amazon- ian forests. Science 299:241–244

► Watson JR, Hays CG, Raimondi PT, Mitarai S and others (2011a) Currents connecting communities: nearshore community similarity and ocean circulation. Ecology 92: 1193–1200

► Watson JR, Siegel DA, Kendall BE, Mitarai S, Rassweiller A, Gaines SD (2011b) Identifying critical regions in small- world marine metapopulations. Proc Natl Acad Sci USA 108:E907–E913

► Watson JR, Kendall BE, Siegel DA, Mitarai S (2012) Chang- ing seascapes, stochastic connectivity, and marine meta- population dynamics. Am Nat 180:99–112

► Whittaker RH (1960) Vegetation of the Siskiyou mountains, Oregon and California. Ecol Monogr 30:279–338