Habitat-Diversity Relations between Sessile Macrobenthos and Benthic Copepods in the Rocky Shores of a Marine Protected Area

: In rocky shore systems, sessile macrobenthic assemblages may act as “ecosystem engineers” for many smaller benthic organisms. Thus, the inﬂuence of macrobenthic coverage on the diversity and assemblage structure of the harpacticoid copepod fauna was investigated in the rocky shores of a Marine Protect Area (MPA) in the Ligurian Sea (NW, Mediterranean Sea). Two sampling sites were investigated in two seasons at three different depths on both sub-vertical and inclined reefs. A total of 61 species of copepods mainly represented by Miraciidae, Laophontidae, Longipediidae and Thalestridae were found. The complex micro-topography of these substrata provided a wide variety of niches for many species with different lifestyles that suggests the important role of rocky shores to ensure the functioning of coastal ecosystems. The harpacticoid assemblage structure seemed mainly inﬂuenced by season and depth. The temporal spread observed is likely one of the underlying mechanisms of niche segregation that allows many species to co-occur in this speciﬁc environment along with a subordinate spatial segregation corresponding to the depth gradient. The results seem to support the hypothesis that the different species composition of the “ecosystem engineer” (and consequently its structure changes) are relevant in structuring the copepod assemblages. The comparison with previous data on general meiofauna underlines that higher surrogacy of the taxonomic identiﬁcation could be used to study rocky shore communities, but the rich diversity that these systems host can only be understood at the lower taxonomic levels. The same holds for future evaluations of impact of environmental changes (including MPA regulations) on meiofaunal assemblages. 2-way crossed SIMPER (Similarity Percentages) carried out to test depth × period interactions. The test shows the percentage dissimilarity between copepod assemblages of the factor


Introduction
Rocky shores host a heterogeneous array of large sessile organisms (e.g., algae, sponges, cnidarians, bryozoans and tunicates) and occurrence of biogenic material (i.e., skeletons of living and dead organisms, shells, wood and rocky clasts) that create a high structural complexity positively related to the biodiversity and, therefore, assumes relevant implications in the conservation of marine coastal systems [1][2][3]. The concept of "ecosystem engineer" was introduced by Jones et al. (1994,1997) [4,5] and defined as a group of organisms that may create, modify or maintain habitats (or micro-habitats) by producing physical state changes in biotic and abiotic variables which, either directly or indirectly, control the food availability for other species. Therefore, rocky-shore-forming

Study Area and Field Sampling
The study was conducted in the Portofino MPA (Ligurian Sea, NW Mediterranean Sea), a marine protected area established in 1999. It is included in the European Natura 2000 Network and is a site of Community Importance [30]. The two sampling sites (Paraggi, P and Aurora, A), both located within the partial reserve zone (Zone C) of the MPA, were selected for the overall low level of human disturbance (Figure 1). In each site, two slopes were selected: sub-vertical (>70°) (1) and inclined (30° to 60°) (2). For each site, sampling stations were located at 5 m (a) and 10 m (b) depth on both sub-vertical and inclined slopes, and at 20 m (c) depth only on inclined slopes (the full list of stations was reported in Table 1). The unbalanced design was due to the geomorphological constraints of the study sites, where sub-vertical walls were absent at depths below 15 m. The sampling activities were carried out during the summer 2012 (i.e., July 2012) and winter 2013 (i.e., January 2013). Detailed information about the sampling routine is reported in Losi et al. (2018) [9]. During the sampling, the composition of the macrobenthic cover was assessed through visual census (quadrat square frame: 50 × 50 cm, divided into 25 sub-squares of 10 × 10 cm each) identifying taxonomically only species that are easily recognizable underwater [31,32].
Meiofauna was collected using a suction sampler (air-lift) designed by Bianchi et al. (2004) [33] to sample hard bottoms. The sampler was equipped with a 38 μm mesh filter, adequate to trap small metazoans (meiofauna). The sampled surface was delimited by a 20 × 20 cm square frame and in order to collect the benthic meiofauna on a standardized area. The spout was repeatedly passed over the rocky substrate, sucking up the surface area delimited by the frame. In each vertical and inclined slope of the two sites (Aurora and Paraggi), three replicates were collected at each depth (i.e., 5 and 10 m in the vertical Detailed information about the sampling routine is reported in Losi et al. (2018) [9]. During the sampling, the composition of the macrobenthic cover was assessed through visual census (quadrat square frame: 50 × 50 cm, divided into 25 sub-squares of 10 × 10 cm each) identifying taxonomically only species that are easily recognizable underwater [31,32].
Meiofauna was collected using a suction sampler (air-lift) designed by Bianchi et al. (2004) [33] to sample hard bottoms. The sampler was equipped with a 38 µm mesh filter, adequate to trap small metazoans (meiofauna). The sampled surface was delimited by a 20 × 20 cm square frame and in order to collect the benthic meiofauna on a standardized area. The spout was repeatedly passed over the rocky substrate, sucking up the surface area delimited by the frame. In each vertical and inclined slope of the two sites (Aurora and Paraggi), three replicates were collected at each depth (i.e., 5 and 10 m in the vertical wall and between 5 and 20 m in the inclined one, see Table 1) after the visual census of the sessile macrobenthos over the same sampling surface.

Meiofaunal and Harpacticoid Assemblages
After the collection, meiofaunal samples were treated with 7% MgCl 2 to promote relaxation of detachment of specimens from surfaces followed by a fixation with 4% buffered formaldehyde, staining with Bengal Rose (0.5 g L −1 ) and storage in 50 mL tubes until laboratory analyses [34]. Once in laboratory, all tubes were sonicated to better detach meiofauna from macroalgae (TRANSONIC LABOR 2000, 3 times for 1 min with 30 s intervals). Then, the meiofaunal fraction was obtained by sieving through 38-1000 µm mesh sieves and centrifugation using a silica gel gradient (Ludox HS40, 1.18 g cm −3 ) [35]. Following De Troch et al. (2001) [36], one hundred individuals or the entire assemblage of harpacticoids (if number was less than 100) were picked out at random from each sample using a stereomicroscope (Leica G26). These obtained specimens were stored in a 75% ethanol and 1% glycerine solution that allows both the maintaining of the body flexibility and prevention of the possible loss of important taxonomical characters such as setae and bristles. Then, each specimen was mounted on a microscopic slide in 100% glycerine, with a small glass support (non-permanent slides) to allow turning of the specimen when observed in toto under a light microscope equipped with Nomarski optics (Optiphoto-2 Nikon) [37,38]. Some specimens were dissected in order to better observe morphological characters otherwise poorly visible [39]. All copepod harpacticoid adults were identified using the following identification taxonomic keys: Lang (1948Lang ( , 1965 [40,41], Huys and Boxshall (1991) [39], Huys et al. (1996) [42], Boxshall and Hasley (2004) [43]. The specimens were identified to species level, when possible, or as morphospecies (indicated as sp.1, sp.2, sp.3, etc.).
Possible synonymies and species distributions were checked using the online site World Register of Marine Species (WoRMS) [44], and the check-list of copepod harpacticoids published by Todaro and Ceccherelli (2010) [45] for the distribution in the Italian seas. Species diversity was reported in terms of richness (number of species), Shannon diversity (H') and Pielou equitability (J', both log2).

Statistical Analyses
Effects of the factors Period (Summer and Winter), Site (Paraggi and Aurora) and Depth (5 m, 10 m and 20 m) on copepod assemblage composition were explored using Water 2021, 13, 1020 5 of 15 several multivariate analyses: the significance of the differences was tested by means of a 2-ways Analysis of Similarities (ANOSIM, Bray-Curtis dissimilarity, data were square-root transformed) considering the following crossed factors: site × period; depth × period and depth nested in sites. The Similarity Percentages-species contributions (SIMPER, 90% cut-off) routine, based on the Bray-Curtis dissimilarity matrix, was performed to measure the level of dissimilarities between factors (i.e., 2-way crossed SIMPER, factor interactions: site × period and depth × period) and find which species contributed mostly to the observed dissimilarities. A non-Multi-Dimensional Scaling (nMDS) performed with a Bray-Curtis Similarity (data square root transformed) was finally used in order to visualize the differences between factors. All statistical tests were performed using PRIMER 6 software package [46].
Principal component analysis (PCA) was carried out to investigate the relationships between copepod fauna and sessile macrobenthic covering. In order to simplify the interpretation of the PCA diagram, the macrobenthic species were grouped at higher taxonomic levels and the various abiotic components (sediment, bare rock, and coarse detritus) in a unique group, namely "Abiotic", while only the harpacticoid species accounting for >2% were considered. All the data were log-transformed prior to the analysis. Data on the sessile macrobenthos (i.e., coverage of sessile macrophytobenthic and macrozoobenthic taxa and abiotic components) were added as primary (or active) variables to the PCA, while data on copepod fauna were entered as supplementary variables without contributing to the overall variance explained by the analysis (see [47] for further details). This last test was applied to provide insights into the possible influence of the sessile macrobenthic cover on the harpacticoid species and univariate measures. The PCA analysis was carried out with the STATISTICA 8 package.
The significance of the univariate variables such as the number of copepod species, Shannon-diversity and Pielou-evenness were checked using the Analysis of Variance (ANOVA). The same statistical design applied for the ANOSIM analysis was applied on ANOVA test. Data that did not meet the normality and homoscedasticity assumptions (Kolmogorov-Smirnov and Levene tests) were log (1 + x) transformed. Instead, Tukey's test was applied for the multiple comparisons when significant differences (p < 0.05) were detected (SPSS software, v.17).

Results
The meiobenthic assemblage was characterized by a high number of taxa (30 taxa), and nematodes and copepods were the most abundant components representing on average 35% and 30% of the whole assemblage, respectively [see 9 for details]. The largest copepod contribution to the meiofaunal assemblage was recorded in summer (50%) in contrast to the low share in winter (17%). The structure of the macrobenthic assemblage was documented in Losi et al. (2018) [9] and reported in Table S1 and exhibited significant variations between the two seasons (p < 0.05), and in particular, the taxa composition between summer and winter differed in both sites (p < 0.05) (see [9] for details).
The complete list of copepod species along with their global and Italian distribution is reported in the Table S2. The copepod assemblage showed a total of 61 putative species belonging to 43 genera and 22 families. Among the species found, Acutiramus brevicaudatus  Table S2).
The two-way crossed ANOSIM carried out assessing site × period interactions revealed significant differences only in relation to the temporal variations (ANOSIM, R 2 = 0.24; p = 0.021), while there were no significant variations of the assemblage structure between Paraggi and Aurora (ANOSIM, p > 0.05). However, when nested ANOSIM was carried out considering the interactions between depth and site, a significant difference between depths was detected (ANOSIM, R 2 = 0.28; p = 0.017), while no significant differences were detected in the interaction depth × period (ANOSIM, R 2 = 0.18; p = 0.056). A visualization of the factor effects (namely sites, periods, depths) on the structure of the copepod assemblage was obtained by means of nMDS (Bray-Curtis similarity, square root transformation Figure 2).
When the site and period factors were analyzed with a crossed SIMPER routine (site × period), they showed a level of average dissimilarity higher between periods (57%) than sites (54%  Table S5).
As documented in the Figure 3, both sites showed very similar values of all univariate measures (total percentage of copepods, number of genera, Pielou and Shannon indices). Additionally, the two periods showed comparable values of richness and diversity indices, while the contribution in terms of abundance of copepods was higher in summer. Instead, a slight increase in number of species and diversity was detected in relation with the depth and, in particular, in the stations located at a depth of 10 m (Figure 3). However, ANOVA did not find significant differences of the univariate measures between the factors analyzed (ANOVA, p > 0.05).   Table S5).
As documented in the Figure 3, both sites showed very similar values of all univariate measures (total percentage of copepods, number of genera, Pielou and Shannon indices). Additionally, the two periods showed comparable values of richness and diversity indices, while the contribution in terms of abundance of copepods was higher in summer. Instead, a slight increase in number of species and diversity was detected in relation with PCA was carried out to analyze the possible influence of the rocky-shore-forming macrobenthos on the copepod distribution. The first two factors of the PCA explained the 49% of the total variance (PC1 27%, PC2 23%, respectively) (Figure 4a). The primary (active) variables of the macrobenthic cover that contributed mainly to the PC1 were Annelida, Ascidiacea, barren substrata, Bryozoa, Rhodophyta and Mollusca, while Ochrophyta, Cnidaria, Porifera, Chlorophyta and turf appeared primarily associated to PC2 (Table 2).  Table 2). The plot of the cases is reported in the Figure 4b. the depth and, in particular, in the stations located at a depth of 10 m (Figure 3). However, ANOVA did not find significant differences of the univariate measures between the factors analyzed (ANOVA, p > 0.05). PCA was carried out to analyze the possible influence of the rocky-shore-forming macrobenthos on the copepod distribution. The first two factors of the PCA explained the 49% of the total variance (PC1 27%, PC2 23%, respectively) ( Figure 4a). The primary (active) variables of the macrobenthic cover that contributed mainly to the PC1 were Annelida, Ascidiacea, barren substrata, Bryozoa, Rhodophyta and Mollusca, while Ochrophyta, Cnidaria, Porifera, Chlorophyta and turf appeared primarily associated to PC2 ( Table

Discussion
Complex ecological interactions between macrobenthic sessile organisms and meiofauna exist in hard substrata assemblages, but the scarce information available on this relationship hampers a comprehensive understanding of their functioning and a possible resilience of these environments [10]. The concept of the complexity of the habitat structure was firstly advanced on bird diversity by Mac Arthur and Mac Arthur (1961) [48] and then it has been applied to aquatic systems. The habitat heterogeneity is clearly a driver for many faunal groups, and ecosystem engineers can greatly amplify it. A large proportion of macrophyte species is recognized as ecosystem engineers' organisms being able to affect the availability of resources to numerous other species by modifying, maintaining, and creating habitats [49,50]. Gee and Warwick (1994) [51] tested whether there was a relation between meiofauna and fractal dimensions in four species of macroalgae and observed that an increasing level of the complexity of algae was positively associated to the meiofaunal biodiversity and led to different assemblage structures. From leaves to stalks of macrophytes, a change of the fauna was documented even at higher taxon level [52]: alga branches were found to be mainly populated by harpacticoids, while the basal part by nematodes [53,54]. Many species of polychaetes (e.g., Sabellariidae) forming large conglomerates of sandy tubes are also recognized as ecosystem engineers and promoters of zoobenthic diversity [55]. Ataide et al. (2014) [56] documented the high significant meiofaunal richness associated to the reefs of Sabellia wilsoni and they supported the theory of Jones et al. (2010) [57] that the settlement of ecosystem engineers in marine bottoms establishes directly or/and indirectly new abiotic conditions; in particular, the high quantity and quality of the trophic resources are certainly reflected as cascading in both benthic and planktonic compartments. Similarly, Passarelli et al. (2012) [11], who mimicked an infrastructure of polychaete tubes, detected that the microphytobenthic biofilm was promoted including the production of the associated extracellular polymeric substances. Thus, it was hypothesized that ecosystem engineers can significantly influence both meiofauna and macrofauna by the creation of complex trophodynamics that raises interesting implications for the maintaining of the ecosystem functioning as well as coastal ecosystem management on the larger scale.
Among meiobenthic organisms, crustaceans were found to be the second most abundant taxon in the MPA of Portofino [9] with a faunal composition comparable to other hard bottoms of the Mediterranean Sea [6,58]. The copepod fauna was dominated by the families Miraciidae, Laophontidae, Longipediidae and Thalestridae that are typically found in phytal habitats [13]. It was interesting to observe that some of the species found in Portofino matched with a species list reported from the Indian Ocean that underlines a high occurrence of cosmopolitan species in the study area and the existence of isocommunities in similar habitats, here hard rock substrates vs. dead coral debris [16].
Copepod species may show numerous and specific adaptations to the habitats (i.e., a different shaped body, anchoring structures), therefore, it is frequently observed not only a change of the faunistic composition, but also in morphological adaptations from the basal to the apical part of the macrophytes [14][15][16][17][18]52,59]. This is well-observable also in copepods from Portofino: Laophonte Philippi, 1840, Esola Edwards, 1891, Porcellidium Claus, 1860 and Paralaophonte Lang, 1948 exhibit legs equipped with claws useful for clinging on epiphytic material and they are, in fact, adapted to live as epibionts on macroalgae (see the dorsoventral flatten habitus of Porcellidium Claus, 1860) especially in environments dominated by strong currents [60]. Thalestris Claus, 1862 species show prehensile maxillipeds and robust hooked appendices to swim and colonize similar typology of algae [14,15,54]. Furthermore, species belonging to the genera Ameira Boeck, 1865 and Amphiascus Sars, 1905, with a clear benthic lifestyle, contributed to a more diverse meiofauna community in the MPA of Portofino.
They occur in the basal part of the macrophytes where sediment accumulates and as such an environment for benthic species is created. Indeed, the species of these genera are generally characterized by an elongated and cylindrical body shape that allows them to better adapt to an interstitial life. However, the occurrence of some eurytopic species, such as Ectinosoma Boeck, 1865 was also revealed [16,52,61]. The presence of both epifaunal and sediment-dwelling taxa is related to the complex micro-topography of the environment, which includes a great variety of niches suitable to host a taxonomical as well as functional diversified meiofaunal assemblage [16,62,63]. Indeed, high stocks and a high diversity of food are present in these types of environments leading to the development of a high variety of selective feeding modes of harpacticoids and the coexistence of species that exploit the same food resources [15].
Assemblage structure revealed more differences in benthic copepods than in general meiofauna [9]: the former, in fact, showed significant changes between both seasons and depths, while meiofauna only temporal variations. In many geographical regions, temporal variations appear to be a primary factor that structures the meiobenthic assemblages of rocky shores. In particular, the biological cycles of epiphytic species (seasonally controlled) seem to be one of the main drivers e.g., [6,[64][65][66] that appears to fit with our PCA results and previous observations on the whole meiofaunal assemblage. Losi et al. (2018) [9], in fact, revealed a significant relation between the cycles of growth and decay of many macrophytes in Portofino substrata. This could mean that the different species' composition (and consequently structure change) of the "ecosystem engineer", seasonally controlled, is the relevant factor influencing the meiobenthic copepods and not temperature variations per se. Overall, the same pattern was found at higher taxon level [9], so we can conclude that higher taxon surrogacy can be applied although the identification at the species level adds extra ecological information [67].
Barren substrata, in the first factor axis of PCA, appeared negatively related to all faunal parameters and copepod species. The macrobenthic species that are more important in PC1 are Annelida belonging to Spirorbidae and Serpulidae families as well as Rhodophyta. The role of polychaeta tubes as ecosystem engineers has already been pointed out above, but also Rhodophyta may attract different crustacean taxa. For instance, Buschmann (1991) [68] documented that crustacean amphipods consume a high amount of Rhodophyta cystocarpic tissues. It is also reported that many copepod species in their different life stages show a very close association with the medullary tissues of both brown and red algae, but Rhodophyta generally contain the most heavy and routine hosts [69]. Thalestridae and Dactylopusiidae are usually recognized as the common families associated to Rhodophyta, but no clear relationship was observed in the present study with these families.
The highest abundance of copepods as well as the larger occurrence of phytal taxa in the study area was found in summer and it was primarily associated to Ochrophyta and secondly to Chlorophyta coverage (Figure 3a). In winter, a general decrease in all harpacticoid species was recorded, likely in relation to the increasing percentage of barren substrates and likely lower food quantity [18]. Furthermore, a negative correlation between copepod diversity and turf percentage was observed unlike that which has been observed for other crustacean species (i.e., ostracods) [10].
Among the most abundant species in the MPA of Portofino, there is Longipedia coronata Claus, 1862 that showed the highest abundances especially in summer. It is known to be a migratory species that is able to live indifferently as epibiont of macroalgae or in the water column [14]. This species has a cosmopolitan distribution and has already been reported for Italian waters (see [70] for review). Another cosmopolitan species that was very abundant in the study area was Laophonte cornuta Philippi, 1840. It was found in several types of substrata from shallow to deep zones and in relation to different types of algae or invertebrates [70].
PCA showed an association between Diarthrodes ponticus (Baird, 1845), barren substrates and Cnidaria, which seems apparently in contrast with the ecological notes on this species reported for soft substrata and phanerogames [71]. However, other representatives of the genus Diarthrodes Thomson, 1883, namely D. nobilis (Baird, 1845), were found in low abundances (<than 2% and so not reported in PCA matrix) in the study site and are known for their ability to produce mucus, an easily degradable substrate for bacteria on which D. nobilis feeds [72][73][74]. This strategy, potentially shared by other congeneric species, might allow it to survive even in barren habitats. The genus Ameira Boeck, 1865 was positively related to Porifera. Indeed, Ameira Boeck, 1865 is known as a sponge dweller [12,15]. The genera Longipedia Claus, 1863, Esola Edwards, 1891, and Phyllothalestris Sars, 1905 were more abundant in the summer period, while Laophonte Philippi, 1840, Porcellidium Claus, 1860, Thalestris Claus, 1862 and Paralaophonte Lang, 1948 were the dominant groups of the phytal assemblage during the winter. This seasonal spread is one of the underlying mechanisms of niche segregation that allows many species to co-occur in this specific environment. Next to a temporal niche segregation, the data also showed a clear spatial segregation of the copepod species along the depth gradient. As mentioned before, a higher taxon surrogacy could be applied to this dataset, but yet this kind of information obtained at the genus level would be missed in case of an analysis restricted to the higher taxa. Therefore, the rich diversity that these environments harbor can only be understood at the lower taxonomic levels. The same holds for future evaluations of impact of environmental changes (including MPA regulations) on meiofauna communities.

Conclusions
The present investigation is one of the few studies on the harpacticoid copepod fauna of rocky shores. A relevant role as "ecosystem engineer" of the sessile macrobenthic assemblages that concurred to maintain a high taxonomical and functional diversity of the harpacticoid species and likely to ensure the resilience of these environments was demonstrated. Season and depth appeared the most important factors influencing the assemblage structure likely because they are involved in mechanisms of niche segregation. It was noticed that the significant temporal variations in the harpacticoid assemblage structure were mainly controlled by different species occurrence and consequent change in the 3D habitat structure of the "ecosystem engineer" and not by temperature variations per se. The comparison with previous data underlines that a higher taxon surrogacy could be used to investigate these coastal systems, but, at the same time, it represents a loss of important information that could help in the evaluations of human impact in MPA and rocky shores.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w13081020/s1, Table S1: Coverage of sessile macrophytobenthic and macrozoobenthic taxa and abiotic components at the investigated stations. Table S2: List of the copepod species found in Portofino Marine Protected Area (MPA) (Ligurian Sea) along with their global and Italian distribution [according to 44 and 33]. Table S3: Percentages of the copepod species detected in the sampling stations of the Portofino Marine Protected Area (MPA). Table S4: Results of the 2-way crossed SIMPER (Similarity Percentages) carried out to test site × period interactions. The test shows the percentage dissimilarity between copepod assemblages of the factor analyzed, as well as the species contributing most to the observed dissimilarity. Cut-off for low contributions 90%. Table  S5. Results of the 2-way crossed SIMPER (Similarity Percentages) carried out to test depth × period interactions. The test shows the percentage dissimilarity between copepod assemblages of the factor analyzed, as well as the species contributing most to the observed dissimilarity. Cut-off for low contributions 90%.