Spatiotemporal Dynamics of Mediterranean Shallow Coastal Fish Communities along a Gradient of Marine Protection

The importance of habitat factors in designing marine reserves and evaluating their performance over time has been regularly documented. Over three biennial sampling periods, we examined the effects of vegetated coverage and habitat diversity (i.e., patchiness) on fish density, community composition, and species-specific patterns along a gradient of protection from harvest in the shallow Spanish southern Mediterranean, including portions of the Tabarca marine reserve. With the exception of two herbivores (Sarpa salpa and Symphodus tinca), vegetated cover did not significantly affect fish densities, while habitat diversity was an influential factor across all three sampling periods. Overall, fish density was more positively associated with more continuous vegetated or unvegetated habitats, and was greatest in areas of highest protection (Tabarca II – Isla Nao site). These patterns were usually observed for four abundant fish species (Boops boops, Chromis chromis, Oblada melanura, and S. salpa). Fish community composition was distinct in the most protected portion of the Tabarca reserve, where it was also most stable. Our findings align with previous investigations of the Tabarca reserve and its surrounding areas, and demonstrate its continued effectiveness in conserving fish biomass and habitat. Together with effective management, marine reserves can facilitate greater species abundance, more stable biological communities, and resilient ecosystems.

In the Mediterranean, multiple investigations have demonstrated the effectiveness of marine reserves and habitat factors for increasing depleted fishery biomass and conserving biodiversity, including their influences on spillover potential and their differential effects on ecologically important species [20,[41][42][43][44][65][66][67][68][69][70][71]. Demersal fishes, in particular, are greatly affected by the quantity and quality of habitat encompassed in Mediterranean protected areas, given their strong site fidelities [41][42][43][44][66][67][68][69][70][71][72]. Recommendations to assess reserve-protective effects with those from habitat through multiscale spatial-integrated approaches have also emerged [73][74][75][76][77], especially for evaluating responses of reef and seagrass-associated fish assemblages-the abundances of which have been affected by overfishing and habitat degradation [17,[20][21][22][23]25,27,[75][76][77]. These spatial factors have been partially considered in previous evaluations of the Tabarca marine reserve in southeast Spain, where some positive responses by fish populations and marine communities have been observed, and spillover gradients developed in relation to the degree of protection and habitat distributions [41][42][43][44][67][68][69]72]. Many of these investigations have taken place within deeper or limited portions of the reserve, with only partial consideration of habitat patchiness and vegetated coverage in shallower locations more frequently affected by human use [41][42][43][44][67][68][69]72]. They have also been mostly conducted within and immediately surrounding the reserve during single-snapshot surveys, or with trends examined over short timeframes. Thus, continued studies in shallower, more regularly accessible portions of the Tabarca reserve and its nearby surrounding areas are warranted for additional evaluation of its effectiveness. Additionally, while previous studies have highlighted depletions of commercially important species in southeast Spain, and examined the effects of human pressures on its coastal habitats, their differential effects on protected and unprotected shallow, coastal fish communities have remained less frequently documented [41][42][43][44][67][68][69]72].
Recurring samplings of marine reserves and their surrounding areas are useful to account for changes in protected and nearby unprotected populations, in addition to habitat dynamics [54,57,66,70,78,79]. Long-term monitoring in the Scandola marine reserve (Corsica, northwest Mediterranean) suggests differential stabilization of fish assemblages in more continuous seagrass versus patchy rocky reef habitats [66,78]. These factors also influence community composition throughout protected and unprotected areas, which likewise contribute to reserve stability and long-term trends in population performance [12][13][14][58][59][60]66,70,78]. Therefore, the variable influences of such spatiotemporal factors on reserve functioning warrant continued investigation.
We examined the effects of habitat coverage and diversity on fish community composition and species-specific patterns within shallow areas along a gradient of marine protection. We chose to focus on shallow, coastal fish communities in accessible areas that would potentially be subject to more frequent human use and probable disturbance. Our specific hypotheses were that greater degrees of protection would result in significantly higher densities of fish species, and altered community composition among sites. In addition, we proposed that greater seagrass coverage and lower habitat diversity in a given sampling area (i.e., more continuous habitats) would most likely occur in protected areas, which would favor higher fish densities, as observed in long-term studies of other Mediterranean reserves [66,78].
Approximately 4 km offshore from SP, the Tabarca marine reserve (~1400 ha., established 1986) surrounds Nueva Tabarca Island-of which, the marine community composition and the effects of its protected status have previously been investigated [41][42][43][44]67,68,72,84,85,93,94]. The reserve is divided into three zones with a gradient of protection levels and allowable uses [41][42][43][44]84,85]. These include a deeper, easternmost 100 ha. integral reserve area (Tabarca I; not sampled in this study), where all activities except scientific study are prohibited; a 630 ha. buffer zone containing Isla Nao (Tabarca II, T2; sampled in this study), where selective artisanal fishing gear for pelagic species is allowed; and the peripheral/transitional zone containing Nueva Tabarca Island (Tabarca III, T3; sampled in this study), where selective fishing gear and recreational activities (swimming, diving, boating) are permitted. Further details regarding zoning, management, administration, participation, monitoring, and scientific oversight have been previously published [72,81,84,85,90]. The entire reserve serves as a no-take zone for species targeted by trammel nets, which are banned throughout this protected area [41].

Methods
Fish density and community composition within the Tabarca protected area (T2, T3) and surrounding locations outside the reserve (CH, SP) were investigated using daytime visual surveys (0800-1600) during three separate sampling periods. Initial surveys occurred from 3-10 September 2009 and were followed by subsequent returns to the study sites during 30 August to 8 September 2011 and 2-11 September 2013. These periods were consistently selected for optimal sampling potential of fish communities during end of summer and early fall seasons. Each year, visual surveys were conducted within two coastal unprotected sites (CH and SP; Figure 1) and two areas within the Tabarca protected area: T2 and the southern coast of the island of Tabarca (T3). Equivalent sampling efforts were conducted during 2009 and 2011 (n = 42 transects; Table 1), while the number of survey transects was increased in 2013 (n = 61) to allow for a more comprehensive evaluation of fish community composition. Within each site, initial 2009 survey locations were haphazardly chosen (i.e., by tossing the transect tape measure in a random direction, allowing it to sink, and swimming to it) in depths no greater than 5 m within areas containing a variable mixture of sandy, rocky, or seagrass (Posidonia oceanica dominant) bottoms. Surveys were executed as 50 m-long × 2 m-wide transects, which were swum by pairs of snorkelers along a tape laid out on the bottom in a randomized direction and following a two-minute acclimation period. Transect start and endpoints were additionally recorded by Global Positioning System (GPS) to allow for subsequent return to these same survey locations. Using a dive slate, a given snorkeler recorded the number of demersal and pelagic fish per species observed along the bottom and in the water column within the 100 m 2 transect. Additionally, information regarding the percent composition of bottom substrate (i.e., sand, rock, and species of seagrass or macroalgae) was recorded along transects at every five meters. To ensure consistency across sampling periods, all teams were trained thoroughly and consistently by the same instructors before carrying out visual surveys. Additionally, during trainings, some transects by paired divers were repeated as back-to-back surveys to test performance and accuracy. Total and species-specific fish densities (number of fish m −2 ) were calculated per transect and substrate classification for all three survey years. The fraction of total vegetated (seagrass and macroalgae) cover and a Shannon-Wiener index of bottom habitat diversity accounting for substrate classifications listed in Table S1 were calculated for each transect. To examine the rigor with which we replicated exact transect paths across the three sampling periods, we calculated the coefficients of variation for the fraction of vegetated bottom cover across the three sampling periods for each transect (mean and SD for each transect). Coefficients of variation varied in magnitude, but many were large, indicating significant technological error (GPS inherent offset) when attempting to locate the same transect over the three sampling periods. Given these observations, data were assumed to meet the assumptions for functional independence. Thus, we opted to analyze fish density data for all fishes and the five-most abundant species separately for each of the three sampling periods.
Two-way effects of site (location) with bottom percent vegetated cover and habitat diversity on fish densities were examined through univariate analyses of covariance (ANCOVAs) and post-hoc Newman-Keuls multiple pairwise testing [95,96]. Given the continuous values for vegetated cover and habitat diversity, ANCOVAs allowed for testing their degree of influence on fish densities among locations. Previous studies have shown the utility of using ANCOVAs to examine habitat effects on fish densities and other biological metrics [97,98]. For each ANCOVA model, the homogeneity of regression slopes assumption was examined through the incorporation of an interaction term (i.e., percent vegetated cover*location or habitat diversity*location) into a given model. Among all models (n = 36), no significant interactions were observed except for in one species-specific case. Additionally, spatiotemporal differences in fish community composition were measured using non-metric multivariate dimensional scaling (nMDS), analyses of similarity (ANOSIMs), and similarity percentages (SIMPER) analyses with PAST 3.22 software [99]. Density data were log-transformed to meet assumptions of parametric univariate testing (Levene's test: F 2,142 = 1.291; p = 0.278), and square-root transformed for multivariate analyses.

Results
A total of 40 fish species were observed among sites, with the majority present in P. oceanica or rock bottom habitats at sites T2 and T3, which were within the Tabarca marine reserve (Table S1). For total fish density, no significant effect of percent vegetated cover was observed, while site was significant for all three sampling periods (Table 2). Additionally, a significant effect of bottom habitat diversity on total density was consistently observed during all years (Table 3), while site was also significant in 2009 (p = 0.000) and 2011 (p = 0.021). Post-hoc comparisons showed higher total densities occurring at T2 than at SP (outside the Tabarca marine reserve) during all years when also accounting for percent vegetated cover (Figure 2 and Figure S1), and during 2009 and 2011 when also accounting for bottom habitat diversity ( Figure 3 and Figure S2). When accounting for vegetated cover, total fish density was generally higher at T2 than at CH (outside the Tabarca marine reserve) and T3, except in 2013. During 2011, no significant difference was observed among total densities when only accounting for habitat diversity (Tables 2 and 3, Figures S1 and S2). Similarly, higher densities were also observed in T3 compared to SP during 2011 and 2013 when also accounting for percent vegetated cover. No significant difference in total fish density was observed between CH and T3 during any of the sampling periods.
Some differences across sites were found for abundant species, particularly when accounting for percent vegetated cover (  (Table S1). In 2009, densities for C. chromis were significantly higher in CH and T2 than in SP and T3 when accounting for either percent vegetated cover or bottom habitat diversity. Additionally, a significantly negative influence of habitat diversity on C. chromis density was also observed during 2009 (p = 0.030) and 2011 (p = 0.015; Table 3; Figure 3b). During 2011, C. chromis densities were highest in T2 than all other sites when accounting for percent vegetated cover, while no significant differences among sites were found when accounting for habitat diversity. In 2013, densities were significantly higher in CH than in SP or T3 when accounting for either percent vegetated cover (p = 0.017) or habitat diversity (p = 0.020), and a significant interaction between percent cover and location was also observed (p = 0.023; Table 2).
During this year, a negative relationship between vegetated cover and C. chromis density was observed in CH, while neutral relationships were observed for all other sites (Figure 2a).     In 2011, densities for Boops boops (Bogue) were significantly lower at SP than at all other sites when accounting for percent cover (p = 0.019), with significant effects of habitat diversity observed during 2011 (p = 0.004) and 2013 (p = 0.025; Tables 2 and 3; Figures 2a and 3a, Figures S1 and S2). For all years, B. boops densities were generally highest in Posidonia and rock habitats in T2 and T3 (Table  S1). In 2011 and 2013, Oblada melanura (Saddled Seabream) densities in T2 were significantly higher than all other sites when accounting for percent cover, with significant effects of habitat diversity similarly observed during these years. O. melanura densities were generally highest in T2 and T3 Posidonia habitats (Table S1). Neutral to negative relationships between habitat diversity and B. boops and O. melanura densities were additionally observed (Figure 3a,b). Densities for Symphodus tinca (East Atlantic Peacock Wrasse) in 2013 were significantly highest in T3 among all sites when accounting for either percent cover (p = 0.003) or habitat diversity (p = 0.004; Table 3; Figures S1 and S2). Densities were generally highest in T2 and T3 Posidonia and rock habitats (Table S1). An additional significant and overall positive effect of percent vegetated cover (p = 0.001) was observed for 2009 S. tinca densities (Table 2; Figure 2b). Although no significant differences across sites were observed for Sarpa salpa (Salema Porgy) densities, significant effects of percent vegetated cover (2009, p = 0.014) and habitat diversity (2011, p = 0.027; 2013, p = 0.014) were found. Densities were generally highest in T2 and T3 Posidonia habitats (Table S1). Overall significant negative relationships between habitat diversity and S. salpa densities were observed, in addition to a negative relationship between S. salpa density and vegetated cover during 2009 (Figure 2b). nMDS ordinations (Figure 4) and ANOSIM analyses revealed significant differences in fish community composition among sites within and across years (Table 4). Community composition in T2 consistently differed from that of all other sites across years, while no significant differences in composition were observed between CH and SP or CH and T3. Although no difference in composition was observed between SP and T3 in the period 2009-2011, these sites were found to differ significantly in 2013 (p = 0.0015). Among these sites, differences in the abundance of O. melanura, S. sarpa, C. chromis, S. tinca, and Diplodus annularis (Annular Seabream) contributed most toward composition dissimilarities (Table 5). Additionally, Trachurus trachurus (Atlantic horse mackerel) contributed up to 11.2% of the differential composition observed between T2 and all other sites in 2009, but did not have any contribution during other periods. Differences in composition among SP and T2 in 2013 were additionally influenced (9.5%) by a differential abundance of Thalassoma pavo (Ornate Wrasse). Table 3. Analysis of covariance (ANCOVA) results examining within-year variability related to Shannon-Wiener habitat diversity index and location (site) for density (fish m −2 ) of all fishes and the five most abundant species. Significant values (p < 0.05) and factors for which significance was observed in a given year are indicated in bold. Interactions (Habitat Diversity*Location) between factors that were significant in at least one given year are additionally included below. All other interactions were non-significant (p > 0.05).

Discussion
Examining the influences of vegetated cover and habitat diversity on fish density throughout a gradient of reserve protection illustrated the differential importance of these factors at the species and community level. Although vegetated cover did not significantly affect fish densities apart from the two herbivores (S. sarpa, S. tinca) during the 2009 sampling period, habitat diversity was a commonly influential factor across all three sampling periods. Overall, and as expected, fish density was positively associated with less diverse and therefore more continuous vegetated or unvegetated habitats among sampling areas, and was greatest in areas of highest protection (T2). Among all sites, fish community composition was distinct in the most protected portion of the Tabarca reserve (T2), where it was also most stable. These findings generally support our initial hypotheses, and demonstrate the effectiveness of protection and habitat dynamics toward reserve functioning in shallow coastal areas. They additionally align with previous investigations conducted throughout other portions of the Tabarca reserve and its surrounding areas [41][42][43][44]67,68,72,85].
In our study, negative relationships between fish density and habitat diversity were found for four abundant species (B. boops, C. chromis, O. melanura, S. sarpa) during at least one sampling period. Additionally, these species were often observed at the greatest densities in the most protected areas of the reserve. The importance of continuous habitats and greater protection in promoting the productivity of these species within the Tabarca reserve was previously observed by Ojeda-Martinez et al. [72]. They attributed these results to the faster growth and higher fecundity associated with these small fish species, in addition to more suitable habitats. Additional studies have also demonstrated the importance of complex and continuous habitats in promoting the abundance of these reef and seagrass associated species [42,[102][103][104]. As observed in our study, Symphodus spp. did not benefit as much from high protection, as potentially related to increased predation pressure in the reserve and the importance of rocky substrate and vegetated cover for these species [72,105]. However, other studies have shown seasonal variation in their feeding behaviors and habitat affinities [106,107], suggesting that there may be differential reliance on habitats within and outside the reserve boundaries throughout a given year. Our findings suggest continued effectiveness of the Tabarca reserve protections in fostering more continuous and less fragmented or anthropogenically disturbed habitats, which in turn has allowed for ongoing enhancement of these populations [41][42][43][44]66]. Given their importance to piscivores, increased abundance of Symphodus spp. may also partially facilitate higher densities of mid-and upper trophic level species (e.g., groupers) in the reserve as noted in previous studies [41][42][43][44]67,85].
Our study highlights the ecological consequences of continued exploitation and human activity throughout southeastern Spanish coastal regions. Low fish densities and high patchiness of habitats continue to be observed in unprotected areas as a result of continued fishing pressure, destructive practices, nutrient loading, and desalination [85,[87][88][89][90][91][92]104,105,[108][109][110]. While 22 reserve networks occur in portions of the Spanish Mediterranean, they only make up~8% of the exclusive economic zone [111,112]. Although management actions are working to address some stressors, their effectiveness has been limited to specific locations [85,[87][88][89][90][91][92]104,105,[108][109][110]. In addition, as other ocean uses and the human population continue to increase throughout the Mediterranean, their effects are likely to intensify within shallow coastal habitats, reinforcing the importance of ongoing protections and broader management approaches [62][63][64][65]85,[87][88][89][90][91][92]104,105,[108][109][110]. Given emergent effects of biological invasions and climate-driven changes to Mediterranean ecosystems, ongoing efforts that include spatial protections are urgently needed to preserve the integrity of their habitats, marine populations, and coastal communities. As recently documented, it is likely that in the time since our study, many factors have continued to affect both unprotected and protected regions of the western Mediterranean [62][63][64][65]80,85,[87][88][89][90][91][92][93]103,105,[108][109][110][113][114][115]. While our study was initiated ten years prior to its publication, the information contained in this work provides valuable baselines for which to compare ongoing monitoring of shallow fish communities and their habitats within and surrounding the reserve. We additionally reinforce the observed stability of fish communities and recurring habitat continuity within its most protected shallow areas for which continued investigations are warranted.
Supplementary Materials: Available online at http://www.mdpi.com/2073-4441/12/6/1537/s1. Table S1. Average annual densities (fish m -2 ; ±SE) of all fish species per survey site and habitat classification. Figure S1. Least-squares mean density (fish m -2 ; ±SE) of all fishes and the five most abundant species per survey site. Uppercase letters indicate significant differences (p < 0.05, analysis of covariance examining site and %vegetated cover, with post-hoc Newman-Keuls tests) among sites for all fishes. Lowercase letters indicate significant differences (p < 0.05, analysis of variance examining site and %vegetated cover, with post-hoc Newman-Keuls tests) among sites for a given species. Figure S2. Least-squares mean density (fish m -2 ; ±SE) of all fishes and the five most abundant species per survey site. Uppercase letters indicate significant differences (p < 0.05, analysis of covariance examining site and habitat diversity, with post-hoc Newman-Keuls tests) among sites for all fishes. Lowercase letters indicate significant differences (p < 0.05, analysis of variance examining site and habitat diversity, with post-hoc Newman-Keuls tests) among sites for a given species.