The continental shelf seascape: a network of species and habitats

The diversity of benthic communities on continental shelves is tightly linked to the diversity of habitats. Therefore, considering seascape habitat composition can help to gain insights into the spatial variability of benthic communities and move away from single-habitats approaches. This perspective needs different analytical methods, such as network analysis that enable the study of complex ecological interactions. This work explores the relationships between habitat and benthic species diversity in the Menorca Channel (the Balearic Islands, western Mediterranean). The seascape in the study area is a mosaic of alternating biogenic and sandy habitats that increases the total benthic species richness. Of the 442 benthic species included in the analyses, 286 species are shared by the six habitats identified, contributing to ecological connectivity across the seascape; 73 generalist species inhabit all six habitats simultaneously, however, 156 species are specialists and are linked to a single habitat, particularly to biogenic habitats, which increases specialization and the vulnerability of the species to habitat fragmentation. The network approach shows a tight link between epibenthic species diversity and the distribution of habitats over the continental shelf, providing essential information for optimal conservation strategies that move from a focus on protecting the most sensitive habitats to marine conservation schemes that encompass a diversity of habitats.


Introduction
Global biodiversity loss and the climate crisis urgently call for conservation and mitigation actions based on a peaceful human interaction with nature (Moranta et al. 2021). However, to advance in this direction, ecological studies are required to generate knowledge at temporal and spatial scales relevant to the conservation of ecosystems (Qiu and Cardinale 2020). The limited spatial dimension of ecological studies is more prevalent in marine than in terrestrial research, due to its complexity and limited access (Raffaelli et al. 2005). This is a paradox as marine habitats cover a vast amount of the Earth's surface, are important contributors to global biodiversity (Davies et al. 2007) and perform several key ecological functions (Snelgrove et al. 2018). The limitation in the spatial extent of marine studies is particularly relevant regarding those ecosystems less accessible for sampling, like deep continental shelf habitats (> 30 m depth). However, progress has been made thanks to large-scale sampling techniques, such as video and remote sensing techniques (Lundsten et al. 2009;Clark et al. 2019;Swanborn et al. 2022).
Soft-sediment habitats that characterise large areas of continental shelves can have small-scale diversity of bottom geomorphology (< 10 m, Hewitt et al. 2005) and biogenic assemblages that can contribute to species richness (Brustolin et al. 2021). Biogenic species aggregations are known to greatly contribute to species' diversity on continental shelves (Hewitt et al. 2008; Barbera et al. 2012;de Juan et al. 2013). However, the complex interconnection between biological communities and the diversity of habitats across spatial scales is currently rarely embedded in the marine ecosystem conservation (Balbar and Metaxas 2019), that often only focuses on the most emblematic habitats (Hewitt et al. 2005). In this context, seascape approaches are crucial to identifying spatially representative priority areas for protection that are characterised by a diversity of habitats ecologically connected over the seabed continuum (Ospina-Alvarez et al. 2020b).
Benthic seascape ecology is rapidly emerging as a key topic within the marine research (Pittman et al. 2021). It has major implications for biodiversity conservation, as it considers the heterogeneous seabed composed of a diversity of interacting components (Pittman et al. 2021;Stuart et al. 2021;Scapin et al. 2022). It is well understood that the diversity and spatial structure of habitat types is a driver of biological assemblages (Swanborn et al. 2022). The spatial composition and structure of habitats link primarily to faunal presence; this includes not only the amount of habitat, cover, or biomass but also the habitat fragmentation (Hewitt et al. 2005;de Juan et al. 2014). Fragmentation transforms the spatial configuration of habitats so that the seascape attributes change, including the reduction of total habitat area and the size of habitat patches, which may impact the diversity of associated fauna (Yeager et al. 2020). Therefore, by considering the effects of seascape composition on the variability of biological assemblages, we gain insights into the vulnerabilities and recovery potentials of species and communities (Lundquist et al. 2010).
The increased awareness of the importance of interactions between ecological components at large scales (e.g., Gonzalez et al. 2017;Erös and Lowe 2019), including marine systems (e.g., Ospina-Alvarez et al. 2020a;Manca et al. 2022), has been paralleled by dramatic improvements in our abilities to describe them. These novel improvements have emerged from the field of graph theory, a branch of discrete mathematics well-known for almost 300 years but that has only recently been incorporated into ecological studies. Graph theory describes the topological properties of a network of interacting elements (Euler 1741) and it is a powerful approach to studying complex ecological interactions (Miranda et al. 2013). These interactions can be represented as graphs or networks, which 1 3 are mathematical abstractions that describe relations between a set of nodes. In ecology, these nodes typically represent ecological entities (e.g., genes, species, populations, or habitats), which are connected by links that represent ecological interactions (e.g., genetic flow, nutrient transport, or population connectivity) (Ospina-Alvarez et al. 2020a). Ecological studies applying network analysis have principally focused on trophic interaction, or plant-pollinators, seed-dispersers, and host-parasite (Manca et al. 2022). The application of network analysis to benthic ecology could provide new insights into the interactions among species and habitats across space (Yang et al. 2018;Manca et al. 2022).
In this work, we explore the relationships between habitat heterogeneity and benthic species diversity at the seascape in the Menorca Channel (the Balearic Islands, western Mediterranean), an area where the continental shelf presents a rich seabed with a dominance of biogenic species. This area harbours a diversity of biogenic and sedimentary habitats that rank among the most important ecotones in the region, including aggregations of red and brown algae that alternate with detritic bottoms and can extend down to 100 m, thanks to the transparency of the water around the Balearic Islands (Ballesteros and Zabala 1993). In the Menorca channel, maërl beds (aggregations of unattached coralline red calcifying algae) and other habitat-forming seaweed species (like Osmundaria volubilis and Peyssonnelia rosa-marina, and brown algae, including the endemic Laminaria rodriguezii (Barbera et al. 2012) are present and in good ecological condition due to hydrodynamic conditions and several fishing-exclusion schemes (Ordines and Massutí 2009). The habitat richness that characterises the Menorca channel provides a natural laboratory to study the interaction between species and habitats in continental shelves.
The diversity of benthic communities is tightly linked to the habitat diversity at large scales and to local environmental characteristics, and these factors are likely to interact (de Juan and Hewitt 2011). Under this assumption, we hypothesized that benthic community diversity is higher in areas with a mosaic of habitats, even when considering low-structured habitats like sand against species-rich biogenic habitats. On the other hand, the species composition across a diversity of habitats conditions the potential ecological connectivity between different habitats, considered as the distance to the nearest population (Thrush et al. 2008). Therefore, the proposed link between habitat heterogeneity and ecological connectivity across seascapes will be key to better understanding processes related to the resilience and recovery of benthic communities after disturbance (de Juan et al. 2014;Gladstone-Gallagher et al. 2019). To explore these hypotheses, we analysed scientific survey data from the Menorca channel that included benthic species inventories across sedimentary and biogenic habitats (Moranta et al. 2014). The existing data sets were explored with network analysis to identify the most important habitats in terms of species richness but also in terms of connectivity to other habitats through shared species.

Study area
The study area is the shallow continental shelf (40 to 100 m depth) between Mallorca and Menorca Islands, known as the Menorca Channel (Fig. 1). The channel has a gentle slope in the first 100 m depth and strong hydrodynamics that attenuates sedimentation (Pinot et al. 2002). The seabed is composed of biogenic sediments, with a predominance of sand and gravel (Druet et al. 2017). Previous studies have identified a diversity of biogenic 1 3 habitats including coralligenous outcrops and maërl beds (Barbera et al. 2012). At the time of the study, the area was under the influence of bottom trawling, except for the area around the underwater communication cables that connect the two islands. However, trawling activities have been excluded from large parts of the Menorca Channel after the designation of the area as a NATURA 2000 site in 2016, which highlighted the presence of maërl and habitat-forming algae beds as important site features.

Data collection
The composition of seabed habitats and their associated biodiversity in the Menorca Channel was studied in three research projects, with fieldwork activities conducted in 2009, 2010 and 2011 (Moranta et al. 2014).
The information provided by epibenthic samples, sediment grabs, side-scan sonar, and video transects allowed the identification of different habitat types characterised by the algal density and the sediment grain size (see details in Barbera et al. 2012) (Fig. 1). The present study focuses on six habitat types: (i) maërl aggregations (with less than 50% of maërl coverage); (ii) maërl beds (> 50% maërl coverage); habitats with the dominance of certain species of algae: (iii) Osmundaria habitats and (iv) Peyssonnelia habitats (areas where Osmundaria volubilis or Peyssonnelia rosa-marina, respectively, dominate over any other biogenic species); (v) sand covered with soft algae (over 30 g/m 2 of algae in the samples); and (vi) bare sand (Table 1).

3
Epibenthic fauna samples were obtained at 146 sites distributed over the 6 habitat types (Table 1). Samples were collected with an epibenthic sledge (2 m wide and a 20 mm codend) towed at 1.5 knots (2.78 km/h) for 5-10 min. All samples were standardized at 1 m 2 . The faunal and algal species retained in each sample were identified to the minimum taxonomic level possible, generally to species level (< 1% identified to genera or higher level). Only strictly benthonic species were retained for the analysis to establish distinctive links between the species and habitats; this implied disregarding demersal fish from the epibenthic species inventories that were previously considered by Barbera et al. (2012). A total of 442 benthic species were included in the analyses.
Grab sediment samples and CTD profiles were obtained to quantify sediment grain size and organic content of sediment, and bottom temperature (Table 1). Side-scan sonar was used to provide information on the nature of the sediments (e.g., texture and consolidation, etc.) and the arrangement of seabed features (e.g., sand waves, rocky outcrops, algae or seagrass beds, anthropogenic features, etc.) (detailed methods in Barbera et al. 2012). The hydrodynamic characteristics of the study area were established with the realistic numerical model DieCAST (Dietrich and Mehra 1998), with the ESEOMED operational system (see detailed methods in Ordines et al. 2011). The bottom current velocity was estimated as the average velocity (cm/s) obtained monthly between 2009 and 2010. Trawling fishing effort was assessed based on Vessel Monitoring System registers for the years 2005-2010 (see detailed methods in Moranta et al. 2014). The number of registers (annual density and annual average) per 3 × 3 km 2 was estimated using ArcGIS v. 9.2 modelling tools.
During the 2011 cruise, which covered 56 of the 146 sites, video transects (ca. 900 m long) were obtained with a video-photo sledge system that independently obtains videos of the bottom and high-resolution images every 2 s with a photographic camera. The analysis of images helped to define habitat characteristics and was approached from two different perspectives: a macroscale analysis from video images to assess habitat fragmentation (number and size of vegetated patches) (Fig. 2a); and a microscale analysis performed through photographs taken each 2 s to assess habitat coverage (Fig. 2b). In this case, image analysis of the surface occupied by each type of substrate was done using an image processing program (Image Tool for Windows v.3.0). This analysis was conducted in one out of every ten photographs, that is, approximately 20 photographs per transect. In each transect, we obtained three parameters: algal cover, number, and size of vegetated patches. The algal  (Moranta et al. 2014). The image on the left: sand with algae; the image on the right: an example of the image analysis of a photograph to calculate the algal cover. In this image, there are two types of small-scale substrate, sand clearings between maërl deposits cover was estimated as the mean value of the percentage of the substrate occupied by algae, calculated from the coverage data of all the photographs treated in each transect. From the images of the videos, the number of patches of vegetation and their average size were estimated. The number of vegetation patches was estimated by counting along the route of the transect the changes in habitat or vegetation, whether it changed from sand to patch with vegetation or if it was different. The size of the patches was calculated as the mean of the size of all vegetation patches in a transect. This size was estimated as a relative proportion (%) = (area of the patch/area of the transect) *100 (Fig. 2b).

Connectivity in species composition across habitats
Species-habitat networks aim at capturing the topological structure linked to the interactions among species and habitats (Yang et al. 2018;Manca et al. 2022); the species are represented by nodes and their weighted co-occurrence across habitat types by links. Both elements define the structure of the underlying network ( Fig. 3). To understand the specieshabitat network(s), the principles of bipartite network analysis were applied (Miranda et al. 2013). A bipartite network is a graph where the vertices (or nodes) can be divided into two disjoint sets such that all edges (or links) connect a vertex in one set (e.g., habitat) to a vertex in another set (e.g., species) (Fig. 3). There are no links between nodes in the disjoint sets (i.e., connections are from habitats to species, and not between habitats or between species). The network analysis applied to the epibenthic species' abundance data allows us to identify the most important habitats in terms of species but also in terms of connectivity to other habitats through shared species.
Considering the composition of habitat patches in the study area, where the different patches alternate creating a mosaic of habitats across the seascape (Fig. 1), the six Fig. 3 Simplified representation of a bipartite network where the species and habitats are the nodes, and the links connect habitats through shared species. The input data for the analysis are the distribution of habitats in the study area ("habitat mapping"); benthic community data in the study area ("species inventories", either as species composition or species relative abundances); then, the habitats and species are coded in nodes and the network is constructed based on the species identified in each node and shared between nodes ("habitat-shared species network", e.g., Fig. 6). Note that the strength of the link between nodes can be weighted by the relative abundance of the species. Finally, a bipartite network can be obtained to illustrate several links to species for each habitat and species shared between habitats (e.g., Fig. 4) 1 3 habitat types were not considered as independent sampling units in the network analysis. Therefore, the full matrix (i.e., 146 sites and 442 species) was used to explore the complete network of species-habitat interactions (i.e., each site is represented by a node). A matrix of interactions was constructed so each species or group of species was explicitly linked to the habitat or habitats in which it occurred. The resulting matrix was then represented as a multiple bipartite weighted network denoting the six habitat types and the link between all species (i.e., the link representing the presence and abundance of the species in the different habitats). The resulting weighted bipartite graph is multiple because an element of one set (species) can be connected to more than one of the elements of the other set (habitats) (Fig. 3). A subset of 107 sites had an estimated level of trawl fishing effort (ranked from no-fishing to high effort); this information was used to create a network of species-habitats labelled by the fishing effort rank to explore the impact of fishing activities on the network configuration.
Several parameters can be calculated to characterise a bipartite network, for example, betweenness, strength or modularity (Kivela et al. 2014;Boccaletti et al. 2014), some of which are node or edge centred (i.e., centrality measures) while others provide global information (e.g., nestedness or modularity).
The centrality measures inform on connectivity properties within the network. Higher values in these centrality measures indicate greater 'connectivity'. In our study context, connectivity is represented through the species shared amongst habitats: species and habitats are the nodes, and the links connect habitats through shared species (Fig. 3). The node strength or weighted degree is the number of species per square meter in the habitat type. The edge-strength is the weight of the flow (of species) between two nodes (see detailed definition of centrality metrics in Ospina-Alvarez et al. 2020a). Node-Strength and Edge-Strength were calculated to highlight the relative importance of each of the species (442), sites (146) and habitat types (6), as interconnected nodes (5886 links). The Page's Rank centrality measure was calculated to identify habitat types with high species abundances, and which are connected (through shared species) to other habitat types with also high species abundances. The centrality measures obtained for each habitat were standardised by the number of sites characterised by that habitat (Table 1). The resulting value is a habitat prioritization index that can inform and guide the protection and preservation of sensitive habitats or groups of habitats.
Global network metrics like nestedness and modularity provide information on the network. Nestedness indicates higher or lower species' generality; a nested network includes a truly diverse set of generalists' species that interact with the more specialised ones (i.e., specialists interact with generalists). Modularity indicates species and habitats that are more closely related to each other than to other species and habitats (Marini et al. 2019). For this purpose, an algorithm based on Beckett's modularity measure, a stochastic optimisation technique that identifies modules in a graph, was used. The Beckett modularity-based algorithms have been used extensively in a wide range of ecological studies for the identification of significant levels of species or community aggregation. These measures were chosen because they are related to network stability (Bascompte et al. 2003) and reflect the structure of the interaction between species and habitats.

Effects of biophysical variables on benthic communities
Biodiversity in the study area is potentially linked to the diversity of habitats and other environmental variables: habitat type, the algal and maërl coverage, the number of habitat patches and patch size for each habitat, depth, organic matter content, grain size, bottom current velocity, and fishing effort. To estimate the effect of environmental variables on the variability exhibited by the community metrics, a sensitivity analysis using generalised additive models (GAMs) was applied to epibenthic community data: abundance and biomass of epifaunal species, and the faunal and algal species richness. GAMs allow the identification of diverse types of relationships between biological and environmental variables, including possible thresholds using thin-plate regression splines, they provide confidence intervals for the regression lines, as well as allow visual inspection of the significance of the relationships. The sensitivity analysis consisted of the following procedure: First, we fitted a GAM. Second, for the GAM, a multi-factor analysis of deviance was performed to determine significance using the R statistical software (R v.4.2.2) and its function "anova.gam" and the Chi-squared test. The stepwise approach for the selection of the best model was: (a) running a full GAM including all variables and their interactions with each of the following distribution errors: Gaussian, Poisson and quasi-Poisson; (b) checking whether the frequency distribution of residuals resembles a normal distribution; (c) checking for over-dispersion (the occurrence of more variance in the data than predicted by a statistical model) and autocorrelation; (d) visual examination of scatterplots of residuals versus predicted values (both in terms of the slope of the relationship and in the dispersion of the values); (e) removal of all non-significant (using AIC with alpha = 0.01) single factors and interactions; (f) calculation of the deviance and residual deviance in order to assess the overall effectiveness of the model and the effects of each independent variable on the dependent variable.
The GAMs were applied to the 56 samples obtained in 2011 that were associated with small-scale habitat structures assessed from the video analysis (Fig. 2). With this data set, benthic community metrics could be linked to habitat variability at the small scale (< 1-10 m): habitat type, the algal and maërl coverage, the number of habitat patches and patch size for each habitat (these data were obtained from the video analysis). The package "mgcv" (Wood 2011) in R was used to optimise the amount of cubic spline smoothing. For the continuous variables, the smoothing function based on the cubic regression spline was used.

The network of species and habitats
The network analysis including the species' relative abundance in the 146 sites characterised by six habitat types shows that Peyssonnelia and sand habitats share numerous species with the other habitats (i.e., they have many connections with species also connected to other habitats) (Fig. 4). Peyssonnelia habitat is consistently ranked highest in the centrality measures obtained (Table 2), however, it is necessary to consider the small number of sampling sites within this habitat (n = 9, equivalent to 6.2% of the 146 sites studied, while 1 3 maërl bottoms represent 36.9%.). The centrality measures were prioritised by standardising the measures by the number of sites per habitat (Table 3). After standardisation, we observed that the highest centrality measures were held by maërl and sand-dominated sites, implying these sites have many species shared with the other habitats.
The more central positions within the network are occupied by the nodes with the highest number of links to other nodes (Fig. 5); hence, habitats or species with a lower number of links are found in the outer regions of the network graph. The habitat with the  highest number of shared species with any other habitat was maërl and the two habitats sharing a higher number of species were sand and maërl, followed by Osmundaria and maërl (Fig. 5). Accordingly, maërl is in a central position sharing species with most habitats. Exploring the presence-absence matrix of species shared between habitats, we observe 1 3 286 epibenthic species shared by the six habitats identified; 73 generalist species inhabit all six habitats simultaneously; and 156 specialist species that are linked to a single habitat, particularly to the vegetated habitats (Osmundaria and Peyssonnelia). Accordingly, the global metrics evidence that the network of species and habitats in the Menorca Channel is more nested than expected by chance (nestedness p-value > 0.001), suggesting that the species-rich habitats host common species also found in the less-species-rich habitats. On the other hand, the network is more modular than expected by chance (modularity p-value < 0.001), suggesting that the habitats tend to harbour a unique assemblage of species (Fig. S1). In summary, the species are shared across habitats, although a group of sites share a distinctive group of species that are less common in the other sites: module 2 includes Osmundaria, Peyssonnelia and sand with algae; module 3, maërl and sand (Fig. S2).
This characteristic is also observed in the network graph based on species relative abundances where the distribution of sites characterised by maërl (red circles) and by sandy habitat (orange circles) in the network overlap on the left half of the graph (Fig. 6), suggesting several species are shared by these two habitat types. On the right half of the graph (Fig. 6), there is a group of nodes representing sand-algae sites (yellow circles), Osmundaria (green circles) and Peyssonnelia (pink circles), which share species. There is also the presence of Habitat type Fig. 6 Network graph of species relative abundances across sites (nested within habitats: legend symbols) in the Menorca Channel. The network is shown as a graph where the 442 species (grey points) and 146 sites (coloured points according to habitat) represent nodes, and the edges represent the connections between species and sites. It is possible to find species that connect to multiple sites but not to other species. The size of the nodes represents the weighted degree (strength), based on abundance. The edges' thickness represents a species' abundance at a particular site sand (orange), maërl (red) and maërl bed (blue) sites on the right half of the graph (Fig. 6), but to a lesser extent. The edges, representing the connections between species and habitat (blue lines in Fig. 6), show an abundance-dependent thickness. Very few species show high abundances that are exclusive to a single habitat type (e.g., Ditrupa arietina related to sand, code "DITRARIE"; Aplidium nordmanni related to Peyssonnelia, code "APLINORD"). There is a group of species in the centre (e.g., Pagurus pridauxi, code "PAGUPRID", and Inachus dorsetensis, code "INACDORS") that are common to most sites.
Some Osmundaria and Peyssonnelia sites also emerge as important in the network graph based on Page's Rank (Fig. S3) occupying a central position and indicating that the most influential nodes in the network are characterised by these habitats. In our case study, a high Page Rank represents sites that share a high number of species with other sites that also share a high number of species. The network analysis overlapping the presence of fishing in the bipartite species-habitat network was done to identify emerging patterns related to fishing activities (presence/absence). We observed that the network follows the structure driven by habitat composition and not by trawl fishing effort (Fig. S4).

Drivers of benthic community variability across the seascape
The variance explained by the GAMs was high, above 80% in some cases, with several environmental and habitat metrics contributing to the explained variance (Table 4). The faunal abundance model evidenced that the size of habitat patches had significant positive effects, while the maërl coverage had a negative effect on faunal abundance. The depth, organic matter and the number of patches had significant non-linear effects on faunal abundance. The model for faunal species richness evidenced significant negative effects of the maërl coverage, whereas the habitat patch size and the number of patches had positive effects on species richness. The speed of bottom currents and the average grain size had significant non-linear effects. The model for faunal biomass evidenced significant negative effects of fishing effort, while maërl coverage had significant non-linear effects. The model for algal species richness explained 40% of the variance and the depth and fishing effort had significant negative effects.

Discussion
The diversity of habitats increases benthic species richness in the Menorca channel. From the 442 species identified in the samples, 73 are generalist species found in all the habitats, while 156 are specialists found in a single habitat. Importantly, 35% of the species (i.e., 286 species) were shared by at least two habitats. These shared species potentially allow for the recolonisation of disturbed patches from neighbouring un-disturbed areas (Lundquist et al. 2010;de Juan et al. 2014). In this context, the Osmundaria volubilis dominated habitat emerges as the most vulnerable habitat to fragmentation as it is the habitat with the higher number of specialist species, which limits potential recovery after disturbance. On the other hand, although sandy sites are characterised by low average species richness, the pool of sandy sites shares many species with other habitats, particularly with maërl beds (which are formed by aggregations of free-living red calcareous algae), increasing the ecological connectivity of the area. The standardisation of the data, by considering the weight of the different habitats in the matrix, proves that maërl is the most central (important) structuring habitat in the Menorca Channel. Overall, the combination of "specialist" and "generalist" habitats (Brustolin et al. 2021), like bare sand and areas with Osmundaria in our case study, will increase total species richness at the seascape while maërl beds are acting as the most connected (structuring) habitat. These results are following previous studies that found that diversified habitats, even including unvegetated habitats, harbour higher invertebrate diversity than typically species-rich habitats alone (Barberá-Cebrián et al. 2002;de Juan and Hewitt 2011).
Several statistical methods have been used in the past to analyse and identify benthic communities' patterns, such as multi-dimensional scaling (MDS), which provide an overview of the distribution of similarities between sampling sites. However, these methods are not suitable for classifying species communities from scratch (Thiergart et al. 2014). Modularity is a key structural feature of complex networks. Accordingly, the identification of modules in complex networks has revealed not only the structured organisation of these networks into intricately linked communities, but also the relationship between this organisation and their functioning (e.g. food web structure) and robustness (Pérez-Matus Table 4 Summary of the GAM for epibenthic community metrics in the Menorca Channel Distribution "Poisson" for richness indices, and "quasi-Poisson" for abundance and biomass of fauna. Linear effects were assessed with "z-value" and splines (s(variable) assessed with Chi.sq Significance levels: *** < 0.001; ** < 0.01; * < 0.05; ns, non-significant  et al. 2017). In species-habitat networks, a modular structure can reveal the flow of common species between neighbouring habitats (i.e., equivalent information processing capabilities). This, in turn, can facilitate the understanding of biodiversity loss in certain habitats by observing how other habitats cluster in the same module (i.e., robustness). It is therefore essential to better understand the modularity of species-habitat interaction networks to uncover fundamental patterns in their structure and function. In our case study, the specieshabitat network is "nested" and "modular", implying the species-rich habitat shares many species with other less rich habitats (considering "rich" by the total accumulated species across sites of a habitat type). However, modularity implies that the habitats are characterised by a unique set of species (and their relative abundance). In summary, the benthic community in the study area is characterised by small invertebrates, which are opportunistic and non-habitat-specific. Despite 35% of the species being linked to a single habitat type, the only two species showing a strong connection (edge strength) to a single site are Diatrupa arenaria to a sandy site (these organisms can accumulate in high densities in soft sediments) and Aplidium normandi, an ascidium associated with Peyssonnelia rosa-marina habitats. This network's central position is occupied by hermit crabs and swimming crabs, as those species shared amongst habitats. Biological communities are composed of rare and common species that contribute in different but important ways to ecosystem diversity and function (Ellingsen et al. 2007;Hinz et al. 2021). In ecological networks, this relative contribution is expressed in the 'core-periphery' structure (Miele et al. 2020). A core of generalist, often common, species are central to maintaining network structure because they tend to be associated with many interactions (Kaiser-Bunbury et al. 2010;Miele et al. 2020). These core generalists interact with both specialists and other generalists, and their loss can lead to cascading species loss and changes in the network function (Bascompte et al. 2003;Memmott et al. 2007). In the species-habitat network from the Menorca Channel, the maërl was found to be the habitat with the highest species diversity (i.e., the higher the number of links, the higher the in-degree) and the highest species abundance (i.e., the more often the habitat-species link is, the higher the in-strength). But also, the sites identified as maërl were connected, through the species, with sites with a high diversity of connections (high Page's Rank). Accordingly, the prioritization index based on the centrality measures highlighted maërl as the most important habitat in the species interaction network, followed by sand.
Maërl beds might act as a structuring habitat in the area that increases species connectivity at the seascape, as it is known to be the substrate where other algae like O. volubilis, P. rosa-marina or Laminaria rodriguezii attach and as such, they have been described as "ecological engineers" (Steller et al. 2003;Barbera et al. 2012). Mediterranean maërl beds harbour a surprisingly high degree of trophic groups and species diversity and contribute to enhancing the biodiversity of surrounding habitats (Barbera et al. 2003). Moreover, dead maërl (often called maërl deposits) is predicted to play an important role in the carbon sequestration (Amado-Filho et al. 2012). Both, live and dead maërl, are important nursery areas and feeding grounds for several species of fish and invertebrates (Ordines and Massutí 2009). The "Habitats Directive" of the European Commission mandates the conservation and management of sea-forming species and the SPABIM Protocol of the Barcelona Convention considers Mediterranean maërl beds as eligible for inclusion in national inventories of sites of conservation importance. However, the impact of fishing gear and dredging on the seabed threaten the health of the Mediterranean maërl beds (de Juan et al. 2013). An important milestone in marine conservation in the Mediterranean is a total ban on the use of trawl gear over maërl beds, such as the trawling ban established in a large section 1 3 of the Menorca Channel in 2016, as a result of the LIFE + INDEMARES project (Farriols et al. 2021).
In the Menorca channel, with most of the area at the time of the study subjected to a low to moderate fishing effort compared to other areas of the western Mediterranean basin at the time of the study (de Juan et al. 2009;Barbera et al. 2012), a sub-network including sites with associated fishing effort level (Fig. S4) showed that the species connectivity between habitats was not modified by fishing activities. The use of Generalised Additive Models confirmed these observations as the fishing effort was not the most important driver of benthic abundance and diversity, however, it had significant effects on biomass. This result supports the prediction that moderate levels of trawl fishing remove large fauna whereas the effects on other components of the benthos become more apparent as fishing intensity increases (Jennings et al. 2001;Hinz et al. 2021). On the other hand, in areas with moderate to low fishing intensity, such as the study area, environmental variability and substrate heterogeneity might mask the effects of trawling activities (Barberá et al. 2017). The small-scale habitat heterogeneity had highly significant effects on the biological metrics, with the size and number of habitat patches having positive effects on benthic species abundance and diversity. These findings agree with studies finding that small-scale habitat heterogeneity increases the diversity of soft sediments (Hewitt et al. 2005;Robert et al. 2014), as long as the biogenic habitat develops to a size that allows the maintenance of benthic species populations (de Juan and Hewitt 2011).
Marine biodiversity conservation schemes should move away from partial and key-habitat-focused protection and adopt a broader seascape perspective (Pittman et al. 2021). A seascape approach is particularly relevant for areas with elevated habitat heterogeneity at local and regional levels such as the Mediterranean (Giakoumi et al. 2013). In this context, biodiversity conservation necessarily needs to acknowledge the importance of certain habitats as high diversity areas (e.g., Osmundaria habitat), whereas other habitats play a key role in harbouring generalist species (e.g., sand) or playing a central role in the ecological connectivity (e.g., maërl). Moreover, habitat diversity contributes to the functioning of the area, as certain species are more associated with sandy bottoms (e.g., D. arenaria, a bioturbator) and others with biogenic bottoms (e.g., A. normandy, a suspension feeder that also increases habitat complexity). The novel application of network analysis allows for exploring the drivers of benthic diversity across the seascape in the study area and identifying the habitat composition that maximises the diversity, as essential scientific input to inform optimal conservation strategies.
Novel management perspectives require alternative ecological assessments of the connectivity of the elements of an interaction network. There is a myriad of useful measures for quantifying centrality in complex networks. In population ecology, they can be useful for finding those habitats or species that excel over other habitats or species and increase the connectivity of the interaction network. That is, to highlight the top-k elements: (1) which of the species living in a network of habitats constitutes the species living in the greatest number of habitats and, (2) which of the habitats connected to other habitats by several species constitute the habitats with the greatest number of species. The global network measures (nestedness and modularity) also allow to explore properties of the species such as generalist, specialist and even endemism at regional scales. Additionally, networks analyses offer a convenient framework to explore how different types of disturbance (e.g., environmental change) can propagate through (and possibly be amplified by) ecological interaction links and, more broadly, affect the structure and stability of the biological systems under study. An informative exercise in ecological network analysis is exploring the robustness of networks to incremental species loss, that is to the subsequent removal of nodes (Memmott et al. 2004). The number of nodes to be removed, as well as the sequence of removal, depends on the specific objectives of the simulation. The exploration of the impact of species' removal through the species-habitat network can be a useful approach to assess the vulnerability of an area to external impact and therefore, identify the most effective spatial conservation measures. This work has demonstrated the great potential of network analysis applied to the study of benthic communities at a large spatial scale. Many advances can be made in this field, for example, by incorporating the spatial dimension in the analysis, so that networks can approximate the seascape configuration and its impact on marine diversity. This can provide crucial information for spatial protection schemes that take into account the small-scale heterogeneity of the seascape, and the need for a holistic approach to the conservation of habitats and their species.
Author contributions SdJ, JM and HH conceived the idea; AO, CB and SdJ designed the methodology; CB and JM contributed to the data; AO and SdJ performed the formal analysis; SdJ led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Declarations
Competing interest The authors have no competing interests or other interests that might be perceived to influence the results and/or discussion reported in this paper.

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