Federico Maioli
email: f.maioli@unibo.it.

Benjamin Weigel

Max Lindmark

Chiara Manfredi

Walter Zupa

Isabella Bitetto

Tommaso Russo

Michele Casini

Assessing the overlap between fishing activities and chondrichthyans distribution exposes high-risk areas for bycatch of threatened species

Keywords: Bycatch, VMS, SDM, Adriatic Sea

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Sharks, rays and chimaeras (chondrichthyans) are crucial for maintaining marine ecosystem functioning [START_REF] Myers | Cascading Effects of the Loss of Apex Predatory Sharks from a Coastal Ocean[END_REF]), but at the same time, they are highly vulnerable to fishing due to their slow growth, late maturity and low fecundity [START_REF] Stevens | The effects of fishing on sharks, rays, and chimaeras (chondrichthyans), and the implications for marine ecosystems[END_REF]. Chondrichthyans are primarily caught as accidental bycatch in fisheries targeting more valuable teleost species, although certain chondrichthyan species are periodically targeted [START_REF] Bradai | Elasmobranchs of the Mediterranean and Black sea: Status, ecology and biology[END_REF][START_REF] Ferretti | By-catch of sharks in the Mediterranean Sea: Available mitigation tools[END_REF]. The harmful effects of fishing are further exacerbated by environmental threats, including habitat degradation and climate change [START_REF] Dulvy | Extinction risk and conservation of the world's sharks and rays[END_REF][START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]. These combined factors have led to a widespread decline in chondrichthyan populations over the past few decades [START_REF] Pacoureau | Half a century of global decline in oceanic sharks and rays[END_REF][START_REF] Simpfendorfer | Widespread diversity deficits of coral reef sharks and rays[END_REF] and have placed many species at risk of extinction [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]. Such declines could, in turn, have adverse impacts on marine ecosystems [START_REF] Ferretti | Patterns and ecosystem consequences of shark declines in the ocean: Ecosystem consequences of shark declines[END_REF]. Spatial management of fisheries is nowadays a largely-applied approach to mitigate the adverse impacts of fisheries (FAO, 2023;[START_REF] Mcconnaughey | Choosing best practices for managing impacts of trawl fishing on seabed habitats and biota[END_REF]. In this context, assessing the spatial overlap between chondrichthyans and fishing effort is crucial for understanding the potential impacts of fishing on these vulnerable species, and thus for designing effective conservation and management strategies [START_REF] Jorgensen | Emergent research and priorities for shark and ray conservation[END_REF][START_REF] Queiroz | Ocean-wide tracking of pelagic sharks reveals extent of overlap with longline fishing hotspots[END_REF][START_REF] Queiroz | Global spatial risk assessment of sharks under the footprint of fisheries[END_REF]. However, limited data on incidental catches, landings, and species identification in many regions hinder the understanding of where, when, and which chondrichthyan species and fishing vessels overlap across their ranges, thereby limiting the implementation of targeted management measures [START_REF] Cashion | Official catch data underrepresent shark and ray taxa caught in Mediterranean and Black Sea fisheries[END_REF]FAO, 2022).

In cases where direct data, such as logbooks, offering accurate and reliable information about catches of these species are unavailable, a widely adopted approach involves the use of tracking devices to monitor the movements of fishing fleets in both space and time [START_REF] Bastardie | Detailed mapping of fishing effort and landings by coupling fishing logbooks with satellite-recorded vessel geo-location[END_REF][START_REF] Russo | VMSbase: An R-package for VMS and logbook data management and analysis in fisheries ecology[END_REF]. Specifically, Vessel Monitoring Systems (VMS) have emerged as a valuable tool for tracking and monitoring fishing effort [START_REF] Amoroso | Bottom trawl fishing footprints on the world's continental shelves[END_REF][START_REF] Eigaard | The footprint of bottom trawling in European waters: Distribution, intensity, and seabed integrity[END_REF]. Through satellite-based communication systems, VMS allows for real-time tracking of fishing vessels, thus facilitating the identification of areas with high fishing activity [START_REF] Witt | A Step Towards Seascape Scale Conservation: Using Vessel Monitoring Systems (VMS) to Map Fishing Activity[END_REF]. By integrating VMS data with species distribution information, a more comprehensive understanding of the spatial overlap between fish and trawling activities can be achieved. This integrated approach may help identify areas of high risk for chondrichthyans bycatch and therefore guides targeted conservation efforts. The Adriatic Sea, situated in the central Mediterranean, has a long history of chondrichthyans fisheries and has experienced depletion of several species in the last few decades [START_REF] Dulvy | Extinction vulnerability in marine populations[END_REF][START_REF] Ferretti | Long-term change in a meso-predator community in response to prolonged and heterogeneous human impact[END_REF][START_REF] Jukic-Peladic | Long-term changes in demersal resources of the Adriatic Sea: Comparison between trawl surveys carried out in 1948 and 1998[END_REF][START_REF] Lotze | Historical Changes in Marine Resources, Food-web Structure and Ecosystem Functioning in the Adriatic Sea, Mediterranean[END_REF]. Among the chondrichthyan species commonly caught or observed in this region are the dogfish (Scyliorhinus spp.), brown skate (Raja miraletus), spurdog (Squalus acanthias), thornback skate (Raja clavata), and smooth-hound (Mustelus spp.) [START_REF] Barausse | The role of fisheries and the environment in driving the decline of elasmobranchs in the northern Adriatic Sea[END_REF]Clodia database, 2020;[START_REF] Ferretti | Long-term change in a meso-predator community in response to prolonged and heterogeneous human impact[END_REF][START_REF] Maioli | Influence of ecological traits on spatio-temporal dynamics of an elasmobranch community in a heavily exploited basin[END_REF]. Presently, approximately 70% of the chondrichthyan species in the Adriatic Sea are regionally threatened according to International Union for Conservation of Nature (IUCN) Red List Criteria (www.redlist.org) [START_REF] Soldo | An Annotated Checklist and the Conservation Status of Chondrichthyans in the Adriatic[END_REF]. This basin is recognized as one of the most trawled regions globally [START_REF] Amoroso | Bottom trawl fishing footprints on the world's continental shelves[END_REF][START_REF] Pitcher | Trawl impacts on the relative status of biotic communities of seabed sedimentary habitats in 24 regions worldwide[END_REF], with bottom trawling being a significant and widely practiced fishing activity targeting various demersal fish species (FAO, 2022). The bottom trawling fleet in the Adriatic Sea consists of 1946 vessels, equivalent to approximately 64,900 gross tonnes, with Italy owning around 70% of the total fleet (FAO, 2022). This fishing activity significantly contributes to the incidental capture of chondrichthyans, making it a major concern for the conservation of these species (FAO, 2022). Recognizing the need for action, the General Fisheries Commission for the Mediterranean (GFCM) initiated the MedBycatch project. Its goal is to develop a collaborative approach to improve the understanding of the Mediterranean multi-taxa bycatch of vulnerable species and test mitigation actions (Carpentieri et al., 2021). However, despite these efforts, there is limited knowledge regarding the bycatch of chondrichthyan species, which hampers effective management and conservation efforts in the region (FAO, 2022). In this study, we adopt a novel approach to examine the spatial co-occurrence between chondrichthyan species and commercial fishing practices using the Adriatic Sea as a case study. We integrate data on fishing effort from the Italian fleets obtained through VMS with predicted species distributions derived from geostatistical Species Distribution Models (SDMs) fitted to data from a fishery-independent bottom trawl survey. Additionally, we analyze and map species richness and the proportion of threatened species. This enables us to assess the overlap between this community and bottom trawl fishing activities, providing insights into the potential impact of fishing. Our goal is to identify high risk areas for bycatch, where the potential impact of fishing on chondrichthyan populations and communities is the highest, thus guiding the development of targeted conservation efforts and management strategies to protect chondrichthyan populations in this heavily exploited ecosystem.

Methods

Study area

The Adriatic Sea, located in the Central Mediterranean between the Italian peninsula and the Balkans, is primarily a shallow and eutrophic basin (Appendix S1). It displays distinct morphological variations along both its longitudinal and latitudinal axes. In the northern part of the Adriatic Sea, there is an extremely shallow mean depth of approximately 30 m up to a maximum of 70 m, and the bathymetric gradient along the major axis of the basin is quite weak (Appendix S1). This region experiences a strong river runoff, with the Po River and other northern Italian rivers contributing around 20% of the entire Mediterranean river outflow [START_REF] Hopkins | The structure of Ionian and Levantine Seas[END_REF]. These rivers introduce significant amounts of nutrients, making the Northern Adriatic Sea the most productive area in the Mediterranean [START_REF] Campanelli | The influence of an exceptional Po River flood on the physical and chemical oceanographic properties of the Adriatic Sea[END_REF] and one of the most heavily fished regions in Europe [START_REF] Eigaard | The footprint of bottom trawling in European waters: Distribution, intensity, and seabed integrity[END_REF]. The central zone of the Adriatic Sea reaches depths of about 200m. However, it is also characterized by two depressions known as the Jabuka/Pomo Pits, which have a maximum depth of approximately 270 m. The Southern Adriatic Sea exhibits significant differences between its northern and southern areas. In the northern region, specifically around the Gulf of Manfredonia, a wide continental shelf can be found, with a distance of about 45 nautical miles between the coast and the 200meter isobath. This area is characterized by a smooth slope. On the other hand, in the southern region, the isobaths are much closer together, with a distance of only 8 nautical miles. These distinct morphological variations between the two areas have direct implications for the biocenoses, fisheries resources, and prevalent fishing techniques employed in the region (Anonymous, 2011). The western coast of the Adriatic Sea is generally regular, sandy, and gently sloping, while the eastern coast is characterized by irregularities, numerous islands, and steeply sloping rocky bathymetry [START_REF] Russo | Adriatic Sea hydrography[END_REF]. The water circulation pattern in the Adriatic Sea is cyclonic, with water masses flowing into the basin from the Eastern Mediterranean following a north-westward flow on the east oat and exiting trough a south-eastward flow on the west coast. The circulation and water masses in the Adriatic Sea are strongly influenced by river runoff and atmospheric conditions, which, in turn, impact the salinity and temperature of the water [START_REF] Artegiani | The Adriatic Sea General Circulation. Part I: Air-Sea Interactions and Water Mass Structure[END_REF].

Fishing effort data

Vessel Monitoring Systems (VMS) data consists of a series of consecutive pings (signals) sent by each vessel at regular time intervals. Compared with other tracking devices like the Automatic Identification System (AIS), the temporal frequency of these pings is relatively low, usually with 1-2 hours between successive pings. However, the VMS data offers unparalleled spatial coverage as it relies on the INMARSAT satellite network [START_REF] Russo | Assessing the fishing footprint using data integrated from different tracking devices: Issues and opportunities[END_REF][START_REF] Shepperson | A comparison of VMS and AIS data: The effect of data coverage and vessel position recording frequency on estimates of fishing footprints[END_REF]. The main limitation of VMS is related to its limited coverage of the fleets since VMS is mandatory only for fishing vessel with a length-over-all larger than 12 m, which means that smaller vessels may not be included in the data. For this study, VMS and logbook data for the trawling fleet operating in the Adriatic Sea were provided by the Italian Ministry of Agriculture, Food Sovereignty and Forests (MASAF) within the scientific activities related to the Italian National Program for the Data Collection in the Fisheries Sector (INPDCF). We extracted a subset of VMS data from 538 vessels utilizing bottom otter trawl (OTB) and beam trawl (TBB) gear types, which are two common fishing methods used along the Italian Adriatic coast. These bottom trawlers were identified through cross-analysis involving logbook data and the EU Common Fleet Register (https://webgate.ec.europa.eu/fleet-europa/index_en). This subset of data was then used to assess the fishing activity of these vessels from 2009 to 2021 using the R package VMSbase (version 2.2.1; [START_REF] Russo | VMSbase: An R-package for VMS and logbook data management and analysis in fisheries ecology[END_REF]Russo et al., 2011a;Russo et al., 2011b). VMS pings were partitioned into fishing trips by vessel and interpolated to increase their temporal frequency to 10 minutes (Russo et al., 2011a). The high-frequency interpolated VMS pings were then joined with the NOAA-Etopo1 database [START_REF] Amante | ETOPO1 arc-minute global relief model: Procedures, data sources and analysis[END_REF] through the R package marmap [START_REF] Pante | marmap: A Package for Importing, Plotting and Analyzing Bathymetric and Topographic Data in R[END_REF] to estimate the bottom sea depth for each ping. Finally, to identify fishing activities accurately, fishing set positions were isolated from other vessel states, such as steaming and resting, by applying a combination of speed and depth filters [START_REF] Russo | VMSbase: An R-package for VMS and logbook data management and analysis in fisheries ecology[END_REF]. Thereafter, we filtered the data to include only the period covered by the MEDITS survey (May to September, see below). Subsequently, we computed the yearly aggregated value of trawl fishing effort (i.e. the sum of the fishing hours) for each cell in a 2 x 2 km square grid multiplying the number of fishing set positions by the interpolation frequency (10 minutes). We then filtered this grid by excluding depths not covered by the MEDITS survey, specifically excluding depths < 10 m and > 800 m. Additionaly we excluded the national territorial seas of Slovenia, Croatia, Montenegro and Albania. The shapefiles for the exclusion of the territorial seas was obtained from www.marineregions.org (2023).

Survey data

The complete biological dataset considered for this study consisted of the catches of 5,122 bottom trawl hauls that were performed in the Italian and international waters of the Adriatic Sea during the Mediterranean International Trawl Survey (MEDITS). This survey has been conducted annually since 1994 by the Laboratory of Marine Biology and Fisheries of Fano (Italy) in the Northern and Central Adriatic Sea and by Laboratorio Provinciale di Biologia Marina of Bari (Italy) (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and COISPA Tecnologia & Ricerca (since 2009) in the Southern Adriatic Sea. The survey is generally carried out in the late spring-summer period (May to September, although in some years it has been performed in October-December), and the hauls are located at depths ranging from 10 to 800 m, following a random-stratified sampling scheme based on 5 different depth strata (Anonymous, 2017;[START_REF] Spedicato | The MEDITS trawl survey specifications in an ecosystem approach to fishery management[END_REF]. The GOC-73 experimental bottom trawl is used as the sampling gear, which has a horizontal opening of 16-22 m and a vertical opening of approximately 2.4 m. The codend of the trawl is equipped with a 20 mm side diamond stretched mesh. Further information on the sampling procedures, data collection, and analysis can be found in the MEDITS handbook (Anonymous, 2017). We selected and analyzed a total of 4,216 hauls in this study (Appendix S1), after excluding hauls conducted prior to 1999 and hauls conducted during the autumn period (October to December). The decision to exclude these hauls was driven by two factors. Firstly, the availability of the biogeochemicals variables (see section Environmental covariates and model selection) only reaches back to 1999. Secondly, the coverage of the autumn season has been inconsistent over the years. Ensuring a seasonally homogenous dataset was essential, considering the potential redistribution of chondrichthyan species across seasons [START_REF] Manfredi | Temporal and spatial distribution of some elasmobranchs in the northern and central adriatic sea/distribuzione spazio-temporale di alcuni elasmobranchi in alto e medio adriatico[END_REF]. Moreover, we excluded species scarcely represented in the data (i.e., those that occurred in less than 3% of all the trawl hauls considered). Further, we grouped the common smooth-hound (Mustelus mustelus) and the blackspotted smooth-hound (Mustelus punctulatus) together under the category smooth-hound (Mustelus spp.) because the morphological identification for the Mustelus genus is not straightforward [START_REF] Marino | Resolving the ambiguities in the identification of two smooth-hound sharks (Mustelus mustelus and Mustelus punctulatus) using genetics and morphology[END_REF], and the two species can therefore be misclassified. Additionally, the two species occur together in the Northern Adriatic Sea, display a similar diet [START_REF] Lorenzo | Ontogenetic trophic segregation between two threatened smooth-hound sharks in the Central Mediterranean Sea[END_REF], inhabit similar habitats, and the chance of hybridization is not excluded [START_REF] Marino | New Molecular Tools for the Identification of 2 Endangered Smooth-Hound Sharks, Mustelus mustelus and Mustelus punctulatus[END_REF][START_REF] Marino | Resolving the ambiguities in the identification of two smooth-hound sharks (Mustelus mustelus and Mustelus punctulatus) using genetics and morphology[END_REF]. In total, 10 species/genera were included in further analyses (Table 1).

Species distribution models

We utilized geostatistical generalized linear mixed-effect models (GLMMs) as SDMs to capture the distribution patterns of chondrichthyan species (Table 1). These models combined survey data and environmental covariates, while also considering the influence of spatial and spatiotemporal correlations using Gaussian random fields. Specifically, we modelled the presence/absence of each species using a Bernoulli distribution:

𝑦 𝒔,𝑡 ∼ Bernoulli( 𝑝 𝒔,𝑡 ) (1) logit(𝑝 𝒔,𝑡 ) = 𝛼 𝑡 + 𝑿 𝒔,𝑡 𝜷 + 𝜔 𝒔 + 𝛿 𝒔,𝑡 (2)
𝛼 𝑡 ∼ Normal(0, 𝜎 𝛼 2 ) (3) 𝜔 𝒔 ∼ MVNormal(𝟎, 𝚺 𝜔 ) (4) 𝜹 𝑡=1 ∼ MVNormal(𝟎, 𝚺 𝜖 ) (5) 𝜹 𝑡>1 = 𝜌𝜹 𝑡-1 + √1 -𝜌 2 𝜖 𝑡 , 𝜖 𝑡 ∼ MVNormal(𝟎, 𝚺 𝜖 ) (6)
where 𝑦 𝒔,𝑡 is the binary response variable (presence/absence) for the observation at space 𝒔 and year 𝑡, 𝑝 𝒔,𝑡 denotes the probability of observing the species at space 𝒔 and year 𝑡, 𝛼 𝑡 represents random intercepts by year 𝑡. 𝑿s,t represents the vector of the environmental covariates (see the section Environmental covariates and model selection and Appendix S2) at space s and year t and 𝜷 represents a vector of corresponding coefficients. ω 𝒔 and 𝜖 𝒔,𝑡 are the spatial and spatiotemporal random fields with estimated covariance matrices 𝚺 ω and 𝚺 ϵ , respectively. Spatiotemporal random effects are here assumed to follow a stationary AR1process where ρ represents the correlation between subsequent spatiotemporal random fields.

To overcome convergence issues encountered during the initial analysis, for the deep-sea species (blackmouth catshark, rabbitfish, and velvet belly lanternshark) and the rare species (with a frequency of occurrence ≤ 5%: starry skate, common eagle ray, brown skate, smoothhound) we did not include spatiotemporal random fields (𝜖 𝒔,𝑡 ) in the models and relied exclusively on the spatial random fields (𝜔 𝒔 ).

To optimize computational efficiency, we utilized a predictive process modeling framework [START_REF] Anderson | Black swans in space: Modeling spatiotemporal processes with extremes[END_REF][START_REF] Latimer | Hierarchical models facilitate spatial analysis of large data sets: A case study on invasive plant species in the northeastern United States[END_REF] for fitting all models. In this framework, we employed the Stochastic Partial Differential Equation (SPDE) approximation [START_REF] Lindgren | An Explicit Link between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach[END_REF] and a triangulated mesh (Appendix S3), generated using the R-package R-INLA [START_REF] Rue | Approximate Bayesian Inference for Latent Gaussian models by using Integrated Nested Laplace Approximations[END_REF], to represent spatial random fields. Random effects are estimated at the knots (the vertices of the mesh) and are then bilinearly interpolated to data locations. The optimal knot locations were determined using a k-means clustering algorithm, which minimized the total distance between knots and data points. We selected 200 knots for our models after an initial exploration, as this balance between accuracy and computational time was deemed appropriate (Appendix S3). The models were fitted using TMB [START_REF] Kristensen | TMB: Automatic Differentiation and Laplace Approximation[END_REF] through the Rpackage sdmTMB (version 0.2.1) (Anderson et al., 2022 [preprint]) with maximum marginal likelihood and the Laplace approximation to integrate over random effects. We assessed convergence of the models and ensured that the maximum absolute gradient with respect to all fixed effects was < 0.001, and that the Hessian matrix was positive-definite (Anderson et al., 2022 [preprint]). The model residuals are presented in the Appendix S4-S7.

Environmental covariates and model selection

For the SDMs, we selected environmental covariates based on data availability and their assumed relevance to the distribution of the studied species [START_REF] Maioli | Influence of ecological traits on spatio-temporal dynamics of an elasmobranch community in a heavily exploited basin[END_REF]. We included four covariates: depth (m), seafloor temperature (°C; hereafter temperature), seafloor dissolved oxygen (ml/L; hereafter oxygen), and seabed substrate (categorized into "Sandy mud", "Fine mud", "Sand" and "Muddy sand"). We extracted monthly mean temperature and oxygen values from the Copernicus Marine Service [START_REF] Escudier | Mediterranean Sea Physical Reanalysis (CMEMS MED-Currents, E3R1 system): MEDSEA_MULTIYEAR_PHY_006_004[END_REF][START_REF] Teruzzi | Mediterranean sea biogeochemical reanalysis (CMEMS MED-Biogeochemistry, MedBFM3 system) (version 1)[END_REF], wich have a spatial resolution of approximately 4 km. We obtained seabed substrate and depth data from the EMODnet Seabed Habitats project (https://www.emodnet-seabedhabitats.eu) and the EMODnet Bathymetry project (https://www.emodnet.eu/en/bathymetry), respectively. To integrate these environmental covariates with the trawl hauls data, we associated each trawl haul with the nearest spatial point of the variables. In the case of temperature and oxygen variables, we also matched the trawl hauls data with the nearest available measurements in time. Inter-annual variations were accounted for by including year as random intercept.

We then asked which combination of covariates resulted in the most parsimonious model, as estimated by the Akaike Information Criterion (AIC). To investigate the relationship with depth, for the most common species (spurdog, small-spotted catshark, and thornback skate), we explored two representations to account for bell-shaped relationships: log(depth) + log(depth 2 ), and depth + depth 2 . Additionally, we evaluated whether a linear term alone or a combination of linear and quadratic terms for temperature provided a better fit to the data. Finally, we tested the inclusion of substrate type and of the interaction term between oxygen and temperature.

To address convergence issues encountered during the initial analysis, we narrowed down our model selection process for the deep sea species and for the rare species. For the deep sea species we focused on the functional form of the depth covariate. For the rare species, we focused on the functional forms of both depth and temperature covariates. Specifically, for depth, we considered the following options for both deep sea species and rare species: depth, log(depth), depth + depth 2 , and log(depth) + log(depth 2 ) (Appendix S2).

Moreover, for the small-spotted catshark, the velvet belly lanternshark and the thornback skate, we refitted the best candidate model excluding the random year intercepts, as their effect sizes were small and yielded convergence issues. For the common eagle ray we set Normal(0,10) weakly informative priors for the linear and quadratic terms of the depth covariate, as the initial estimates were unreasonably high, and consequently, we refitted the best candidate model using a penalized likelihood framework.

In total, we compared 16 different models for each common species, 4 for each deep sea species, and 8 for each rare species. To facilitate comparison of effect sizes between the included covariates and to compare them to the marginal standard deviation of spatial variation, we rescaled all covariates to have a mean of 0 and a standard deviation of 1 [START_REF] Schielzeth | Simple means to improve the interpretability of regression coefficients: Interpretation of regression coefficients[END_REF]. Finally, the explanatory power of the models was evaluated by computing the R 2 values [START_REF] Nakagawa | The coefficient of determination R 2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded[END_REF][START_REF] Nakagawa | A general and simple method for obtaining R 2 from generalized linear mixed-effects models[END_REF] and the Area Under the Curve (AUC; [START_REF] Pearce | Evaluating the predictive performance of habitat models developed using logistic regression[END_REF].

Prediction grid

We predicted the probability (p) of chondrichthyan species occurrence onto the same 2 x 2 km grid used for the assessment of the trawl fishing effort, using the best-fitting models. We followed the previously mentioned methodology by matching environmental covariates to the midpoint of each cell. For this particular step, we utilized the average values of temperature and oxygen spanning the core survey period, which extends from May to September. By considering the mean values over this period for each cell, we standardized the effects of temperature and oxygen, making our predictions comparable across different years. We restricted the predictions to the period between 2009 and 2021, as VMS data were unavailable for earlier years.

Species richness

To estimate the species richness, we predicted onto the grid by drawing 1000 simulations from the joint precision matrix of the best-fitting model for each species. This allowed us to obtain an estimate of species richness by summing the number of species in each grid cell [START_REF] Ovaskainen | Joint species distribution modelling: With applications in R[END_REF].

Proportion of threatened species

To characterize the chondrichthyan community as a whole, we calculated the proportion of threatened species. This was done for each cell of the prediction grid, weighing the speciesspecific conservation status by the predicted proportion of each species in each cell. To determine the conservation status of each species, we used the criteria set forth by the IUCN (www.redlist.org), categorizing species as threatened (1) if they were classified as Critically Endangered (CR), Endangered (EN), or Vulnerable (VU), and as nonthreatened (0) if they were categorized as Near Threatened (NT) or Least Concern (LC) (Table 1). Spurdog, common eagle ray and smooth-hound are categorized as threatened species, whereas the remaining species are classified as nonthreatened (Table 1).

Spatial Overlap Index

To evaluate the relationship between the distribution of chondrichthyans and fishing activities, we developed a species-specific Spatial Overlap Index (SOI).

SOI 𝑖,𝑔,𝒔,𝑡 = 𝑝 of presence 𝑖,𝒔,𝑡 × Fishing effort 𝑔,𝒔,𝑡 Max(𝑝 of presence 𝑖,𝑡 × Fishing effort 𝑔,𝑡 ) (7)

The calculation for the SOI involves multiplying the predicted probability (p) of species i being present in a specific grid cell s for a given year t by the corresponding fishing effort (measured in fishing hours) for gear type g during the same year t and grid cell s. The resulting values are subsequently divided by the maximum value obtained from the product of predicted probabilities with fishing effort for gear type g and year t. This rescaling ensures interpretability and constrains the index within a range of [0,1], and puts a stronger emphasis on the spatial overlap patterns. Consequently, we presented the SOI as an average across the most recent years (2018-2021) since we consider this timeframe more relevant from a conservation point of view. Also, since 2018 a Fisheries Restricted Area was established in the Jabuka/Pomo Pits banning demersal fisheries (GFCM, 2021).

Range overlap

We employed a range overlap metric (RO) as defined in [START_REF] Carroll | A review of methods for quantifying spatial predator-prey overlap[END_REF] to quantify the overlap between chondrichthyan species and trawling activities. The RO metric estimates the proportion of a species' area of occupancy where also trawling activities occur:

𝑅𝑂 𝑖,𝑔,𝑡 = 𝐴 𝑖,𝑔,𝑡 𝐴 𝑖,𝑡 ⁄ (8)

Where 𝐴 𝑖,𝑔,𝑡 is the area occupied by the species i and trawled by gear g in year t, while 𝐴 𝑖,𝑡 is the total area occupied by species i in year t [START_REF] Saraux | Spatial Structure and Distribution of Small Pelagic Fish in the Northwestern Mediterranean Sea[END_REF]. Here, occupied cells were defined as cells with ≥ 75th percentile of the species maximum probability of occurrence.

Similarly, to identify cells with high trawling activities, we used a threshold where the fishing effort corresponds to ≥ 75th percentile of the maximum fishing effort for a gear type. Here, we normalized yearly fishing effort within the range [0,1]. This normalization helps address the potential limitation of VMS data as an accurate proxy for overall fishing effort. Due to the pronounced zero-inflation in grid-cell level VMS data, caused by the spatially heterogeneous fishing effort, we calculate percentiles based only on cells with fishing effort > 0. Cells with a fishing effort of 0 are therefore assigned a value of 0. While an alternative approach involves computing percentiles across the entire domain and selecting a higher threshold, we find it more intuitive to interpret the percentiles of presences in this zero-inflated scenario. To ensure the robustness of the chosen threshold, we also explored the sensitivity to different set thresholds.

Species richness and proportion of threatened species overlap with fishing effort

Finally, we explored the chondrichthyans community overlap with fishing effort, using the estimated species richness and the proportion of threatened species. To visualize these relationships, we employed the R-package biscale [START_REF] Prener | biscale: Tools and palettes for bivariate thematic mapping[END_REF] to generate bivariate choropleth maps. To this end, we used fishing effort, species richness, and the the proportion of threatened species as averages across the most recent years (2018-2021).

Results

Fishing effort

Overall, we processed an average of 445 unique VMS identifiers for otter bottom trawlers (OTB) and 93 for beam trawlers (TBB) per year (Appendix S8). Otter trawling was extensively practiced throughout the entire region (Appendix S9). On the other hand, beam trawling was primarily concentrated in the northern regions, extending along the coastline towards the central part of Italy (Appendix S9). Fishing effort (fishing hours in our study) was substantially higher for OTB compared with TBB. The temporal variation over space of OTB was relatively lower (lower coefficient of variation, CV) compared with TBB (Appendix S9), which can be attributed to the lower fleet size of the latter (Appendix S8).

Species distribution models

Considering the entire survey period (1999-2021), the most dominant species were smallspotted catshark and blackmouth catshark being observed in ~ 10% of the hauls, followed by spurdog, velvet belly lanternshark, thornback skate, and rabbitfish being recorded in ~ 8% to 6% of the hauls, while starry skate, common eagle ray, brown skate, and smooth-hound occurred more sporadically (Table 1).

All species models demonstrated a good fit to the data, with conditional R 2 values ranging from 0.54 to 0.93, and AUC values ranging from 0.91 to 0.99 (Appendix S2). Overall, the inclusion of the linear and quadratic terms for depth provided better support to the models (i.e. lower AIC) across species. Specifically, half of the species models showed a better fit using logtransformed depth while the other half performed better without the transformation (Appendix S2). 3 species models exhibited better fitting when using the linear temperature term, another 3 benefited from employing both the linear and quadratic termperature terms, and for 1 species model, the inclusion of temperature was not supported. Among the 3 species for which the interaction term between temperature and oxygen was considered, the inclusion of this term led to better model fitting in 2 cases (Appendix S2).

The models' predictions revealed a distinct spatial pattern in the chondrichthyan community (Appendix S10, S11). Specifically, among the threatened species, spurdog exhibited an elevated average probability of occurrence (p) across the northern and central offshore zones. Similarly, common eagle ray and smooth-hound displayed higher average p values in the northern coastal regions, extending into the offshore area below (Appendix S10).

Conversely, for nonthreatened species, the small spotted catshark, thornback skate, and brown skate predominantly occupied the central offshore areas and the southeastern part of the study area (Appendix S11). Deep-sea species such as the blackmouth catshark, velvet belly lantern shark, and rabbitfish were prevalent in the deeper southern central regions. The starry skate, on the other hand, showed a more prevalent distribution along the central and southern coastal zones (Appendix S11). Species richness exhibited a distinct contrast between the coastal central and southern regions, which were relatively poor in species, and the northern and deeper southern regions, which were rich in species (Fig. 1a). In the northern areas, the relatively high probability of occurrence of spurdog, smooth-hound, and common eagle ray (Appendix S10) contributed to the observed higher species richness. On the other hand, the higher species richness observed in the deeper southern regions was driven by the local high probability of occurrence of deep-sea species such as blackmouth catshark, velvet belly lanternshark, and rabbitfish (Appendix S11). Although the overall level of uncertainty in the mean predicted species richness was relatively low (SD < 0.4), it exhibited higher values in the central and southern offshore areas, suggesting greater temporal variations (Fig. 1b). The proportion of threatened species was found to be highest in the northernmost part of the study area and lowest along the central and southernmost coasts (Fig. 1c). This pattern was influenced by the relatively high local probability of occurrence of spurdog, smooth-hound, and common eagle ray in the northern regions, all of which are classified as threatened species (Appendix S10). The overall level of uncertainty in the mean proportion of threatened species remained relatively low (SD < 0.2), exhibiting lower values in the northern and central coastal regions as well as in the southern and southeastern offshore areas. This suggests that these specific areas experienced lower temporal fluctuations (Fig. 1d).

Spatial Overlap Index

For the threatened species, the analysis of the Spatial Overlap Index (SOI) between species and fishing effort across the most recent years (2018-2021) revealed high mean overlap scores for the spurdog, common eagle ray, and smooth-hound in the northernmost part of the study area, both for OTB and TBB (Fig. 2). Additionally, spurdog exhibited a notable level of overlap with OTB along the offshore regions in the northern-central part of the basin. The SOI for the nonthreatened species is reported in Appendix S12-S13. Despite year-to-year fluctuations, the SOI values remained relatively consistent over time, with closer years showing higher similarity (Appendix S14-S17).

Range overlap

For several species, there was a high potential risk of bycatch extending across a significant portion of their distribution within the study area, as indicated by high range overlap values (RO). Among the species analyzed, starry skate, spurdog, and velvet belly lantern shark were found to share a significant overlap with regions of intensive OTB fishing effort, as indicated by higher RO values (≥ 0.15, averaged across years). Rabbitfish, common eagle ray, blackmouth catshark, smooth-hound, and brown skate exhibited somewhat lower but still significant RO values (< 0.15 and ≥ 0.1, averaged across years) (Fig. 3). On the other hand, common eagle ray, smooth-hound, spurdog and starry skate had the highest overlap with regions affected by TBB activities (≥ 0.1, averaged across years) (Fig. 3). The RO metrics exhibited variations across years, with significant fluctuations observed for certain species, including spurdog, common eagle ray, and smooth-hound, in both OTB and TBB contexts (Fig. 3). Temporal patterns were also evident, with RO values increasing or decreasing together in certain years, likely in response to changes in fishing effort.These findings remained consistent regardless of the threshold used (Appendix S18-S19).

Species richness overlap with fishing effort

In general, the areas exhibiting simultaneous high species richness and high OTB fishing effort were predominantly found in the northernmost section of the basin and in some areas along the outer edge of the continental shelf in the southernmost region (Fig. 4). On the other hand, for TBB, our analysis revealed significant overlap between regions of high species richness and fishing effort in the northernmost area of the basin and the northern central region (Fig. 4). Additionally, we identified two distinct zones in the northern-central part of the basin and in the deeper offshore southernmost areas, characterized by high species richness and comparatively low fishing effort of OTB and TBB (Fig. 4).

Proportion of threatened species overlap with fishing effort

The analysis of OTB revealed distinct patterns in the spatial overlap between fishing effort and the proportion of threatened species. Specifically, in the northernmost regions of the study area, high fishing effort coincided with a higher proportion of threatened species (Fig. 5). However, towards the central and southern portions of the study area, the proportion of threatened species gradually diminishes, in combination with the consistently high levels of OTB fishing activity.

Regarding TBB, a clear spatial pattern emerges. The northern part of the study area exhibits a pronounced overlap between fishing effort and the proportion of threatened species. However, this overlap rapidly diminishes southward (Fig. 5). Moreover, we identified distinct areas in both the northernmost and northern-central parts of the basin, which were characterized by a medium-to-high proportion of threatened species and relatively low fishing effort of OTB and TBB (Fig. 5).

Discussion

Effective marine conservation and fisheries management requires a comprehensive understanding of the spatial overlap between species distribution and fishing activities [START_REF] Douvere | The importance of marine spatial planning in advancing ecosystem-based sea use management[END_REF]. Especially for chondrichthyans, the deficiency of such information is often recognized as a key factor contributing to the failure to implement suitable management measures [START_REF] Cavanagh | Overview of the conservation status of cartilaginous fishes (Chondrichthyans) in the Mediterranean Sea[END_REF]FAO, 2022). Bottom trawl fishing is known to result in incidental bycatch of nontarget species, including chondrichthyans (Carpentieri et al., 2021).

In addition to the direct removal of chondrichthyans, bottom trawl fishing also poses a threat to their habitats. The use of weighted nets dragged along the seabed can cause disturbance and destruction of benthic habitats and removal of the associated biota [START_REF] Ices | Working group on fisheries benthic impact and trade-offs (WGFBIT; outputs from 2022 meeting)[END_REF][START_REF] Sciberras | Response of benthic fauna to experimental bottom fishing: A global metaanalysis[END_REF][START_REF] Watling | Disturbance of the Seabed by Mobile Fishing Gear: A Comparison to Forest Clearcutting[END_REF], which are crucial for the reproduction, sheleter and foraging of many chondrichthyans (e.g. [START_REF] Carrier | Habitat management and closure of a nurse shark breeding and nursery ground[END_REF][START_REF] Heithaus | Habitat use and foraging behavior of tiger sharks (Galeocerdo cuvier) in a seagrass ecosystem[END_REF][START_REF] Simpfendorfer | Utilisation of a tropical bay as a nursery area by sharks of the families Carcharhinidae and Sphyrnidae[END_REF].

In this study, we propose an innovative approach that produces model-based overlap between chondrichthyan distribution and bottom trawl fishing, leveraging high-resolution Vessel Monitoring System (VMS) and the spatial distribution of species derived from geostatistical models fitted to fishery-independent survey data. This information may serve as a critical tool for identifying areas that require protection and implementing conservation measures, such as establishing protected zones or imposing fishing restrictions. Such measures can play a crucial role in safeguarding vulnerable chondrichthyan species and communities, thus promoting their conservation and recovery. We applied this integrated approach to the Adriatic Sea, an highly exploited area [START_REF] Amoroso | Bottom trawl fishing footprints on the world's continental shelves[END_REF] characterized by a long history of chondrichthyan depletion [START_REF] Ferretti | Long-term change in a meso-predator community in response to prolonged and heterogeneous human impact[END_REF] and limited data on their incidental catches (FAO, 2022). We addressed the latter knowledge gap providing novel insights into the potential interaction between chondrichthyans and bottom trawl fishing in this basin.

Our findings reveal significant spatial overlap between chondrichthyans and bottom trawl fishing hotspots in the western part of the Adriatic Sea. This overlap raises concerns about the vulnerability of chondrichthyans, that also includes threatened shark and ray species, to fishing pressure, particularly in regions with high fishing intensity. Our results show distinct spatial patterns of overlap, highlighting both areas of concern and potential refuge for chondrichthyan populations.

We detected two distinct species richness hotspots within the study domain. One hotspot, situated in the northern and central offshore region of the basin, featured primarily threatened species like spurdog, smooth-hound, and common eagle ray, as well as some non-threatened species including thornback skate and small-spotted catshark. The second hotspot, located in the southern and deepest regions, was characterized by the presence of non-threatened deep sea species such as blackmouth catshark, velvet belly lanternshark, and rabbitfish.

In the northernmost part of the Adriatic Sea, we observed a pronounced overlap between the distribution of a high proportion of threatened chondrichthyan species (i.e. spurdog, smoothhound, and common eagle ray) and hotspots of bottom trawl fishing activity. Moreover, for these species, this significant spatial overlap extended over a considerable portion of their habitat, despite some year-to-year fluctuations. These areas of high spatial overlap suggest potential risks for threatened chondrichthyan populations. Conversely, we also identified areas with minimal spatial overlap between high chondrichthyans probability of occurrence and bottom trawl fishing, indicating potential refuge areas. One such area, known as "area dei fondi sporchi" (or "dirty area") by local fishers, located in the northern-central part of the basin, contains relict sand rich in epifaunal organisms that make trawling difficult or even impossible [START_REF] Scaccini | Dati preliminari sulle zoocenosi bentoniche e sulla biomassa in una zona dell'alto e medio Adriatico[END_REF]. Furthermore, our analysis revealed localized differences in spatial overlap between different fishing gears, such as otter bottom trawling (OTB) and beam trawling (TBB), along the central and southern Adriatic coastline. These differences highlight the significance of taking specific fishing gears, and their respective potential impacts, into account when evaluating spatial overlap and formulating effective management strategies. Notably, both OTB and TBB are the focal gears addressed within the Multiannual Management Plan (GFCM, 2019) for demersal species in the Adriatic Sea including hake (Merluccius merluccius), red mullet (Mullus barbatus), common cuttlefish (Sepia oficinalis), Norway lobster (Nephrops norvegicus) and common sole (Solea solea) (Sea Around Us, 2006). This Management Plan combines different measures and among them, fishing effort regimes and spatio-temporal closures. Tailoring conservation measures, within this management context, related to these fishing gears and their interactions with chondrichthyan populations, can optimize conservation efforts and minimize unintended bycatch or habitat damage.

In addition to identifying hotspots of overlap between species or proportion of threatened species with fishing, our model-based approach also allows calculation of range overlaps. This is an important metric along identifying hotspots of overlap, because together with distribution models, it provides information about potential effects of relocating spatial fishing effort. For instance, for species with clear spatial hotspots and low range overlap, relocating fishing effort might reduce mortality from incidental bycatch substantially within the domain. However, for species with high range overlap, simply limiting fishing in "the hottest spot" might have neglible effects, if the fishery over the entire species range and the area with high overlap is large. This illustrates the need to develop a suite of indices before implementing spatial fisheries-restrictions. Understanding where fishing would increase after such an intervention, and how that affects the total mortality rates from fishing, is one of the main challenges in designing protected areas successfully [START_REF] Hilborn | When can marine reserves improve fisheries management?[END_REF][START_REF] Ovando | Assessing the population-level conservation effects of marine protected areas[END_REF].

While our study provides novel insights into the spatial overlap between chondrichthyan distribution and bottom trawl fishing in the western part of the Adriatic Sea, there are certain limitations that should be considered when interpreting the results. Firstly, our analysis relied on the availability and quality of data sources, which may introduce uncertainties and biases. For instance, VMS data is considered a good proxy for fishing effort for the fleet segment larger than 15 meters in overall length, while their reliability decreases for vessels between 12 and 15 meters in length and miss the fishing effort of the fleet below 12 meters. This may result in only a partial representation of fishing activities, potentially leading to an incomplete understanding of the spatial overlap between fishing and chondrichthyan distribution for smaller vessels. Secondly, our data only cover the Italian fleet, while other fishing countries can operate outside the Italian territorial sea. Nonetheless, approximately 70% of the fishing fleet operating in the Adriatic Sea belongs to Italy (FAO, 2022) and therefore we cover the majority of the fleet that should have the major impact on the chondrichthyans in the study area. Additionally, SDMs are subject to inherent limitations, such as uncertainties in species occurrence data, model assumptions, and potential omission of relevant environmental variables. However, the approach described here, calculating overlap on model-predicted probabilities of occurrence also has benefits, primarily in that it is not as sensitive to temporal changes in sampling across in space [START_REF] Thorson | Model-based inference for estimating shifts in species distribution, area occupied and centre of gravity[END_REF]. Lastly, the spatial resolution of our analysis may influence the observed patterns of overlap. While we aimed to use the finest available resolution of the data, the accuracy and precision of VMS data and species occurrence records can vary spatially. Fine-scale variations in both fishing effort and chondrichthyan distribution may not have been fully captured in our analysis, potentially leading to underestimation or overestimation of the true spatial overlap. Despite these mentioned contraints, the approach presented here, to the best of our knowledge, offers the most comprehensive information available regarding the overlap of fishing effort with the chondrichthyan community in the Adriatic Sea.

Our analysis focused on a specific period of the year (from the late spring to the end of summer) due to the lack of consistent autumn and winter fisheries-independent surveys during the last twenty years. Seasonal variations in fishing effort and chondrichthyan distribution may occur though, and therefore caution must be taken in extrapolating our results to other seasons, especially for conservation and management purposes. While spatial overlap between fishing activities and chondrichthyan distribution indicates potential interaction, it does not directly provide evidence of actual biological impacts or quantify the extent of such impacts on chondrichthyan populations. Nevertheless, Mediterranean chondrichthyans are widely acknowledged to be highly vulnerable to trawlers (Carpentieri et al., 2021;[START_REF] Cavanagh | Overview of the conservation status of cartilaginous fishes (Chondrichthyans) in the Mediterranean Sea[END_REF]FAO, 2022). Given the limited availability of direct and quantitative official bycatch data, our results may offer an initial and valid approach to inform conservation efforts and guide precautionary management decisions.

In conclusion, our study introduces an innovative approach that combines VMS data with geostatistical species distribution models, allowing us to analyze the spatial co-occurrence of chondrichthyan species and commercial bottom trawl fishing practices in the Adriatic Sea. Through this analysis, we identified areas of heightened risk for bycatch and thus of potential impacts on chondrichthyan populations, as well as areas that could serve as potential refuge. These findings provide therefore new knowledge that may support the implementation of targeted conservation efforts and the development of effective spatial management measures in the region.

Tables and figures

Figure 1 .

 1 Figure 1. Mean and standard deviation (SD) of species richness (a-b). Mean and SD of proportion of threatened species (c-d). Metrics are computed over the period 2009 to 2021.

Figure 2 .

 2 Figure 2. Mean Spatial Overlap Index (SOI) for threatened species and otter bottom trawling (OTB) and beam trawling (TBB) computed over the period 2018-2021. SOI values of 0 are displayed in gray.

Figure 3 .

 3 Figure 3. Yearly range overlap (RO) between species and otter bottom trawling (OTB) and beam trawling (TBB) fishing gears. A grid cell is considered to be occupied by a species with ≥ 75th percentile of the probability of occurrence, and by trawling activities with ≥ 75th percentile of fishing effort. The year 2017 is omitted due to sampling occurring outside the designated sampling period (May to September; see the section Survey data in Methods).

Figure 4 .

 4 Figure 4. Bivariate plot illustrating the relationship between mean species richness and the percentile of fishing effort for the period 2018-2021 for otter bottom trawling (OTB) and beam trawling (TBB). Dark red areas indicate locations with the highest mean predicted number of species and higest fishing effort. Light red areas correspond to regions characterized by high species richness and low fishing effort. Conversely, dark blue areas represent locations with high fishing effort but a low number of species, whereas gray areas depict regions where both fishing effort and number of species are low.

Figure 5 .

 5 Figure 5. Bivariate plot illustrating the relationship between the mean proportion (prop.) of threatened species and the percentile of fishing effort for the period 2018-2021 for otter bottom trawling (OTB) and beam trawling (TBB). Dark red areas indicate locations with the highest mean predicted proportion of threatened species and fishing effort. Light red areas correspond to regions characterized by a high proportion of threatened species and low fishing effort. Conversely, dark blue areas represent locations with high fishing effort but a low proportion of threatened species, whereas gray areas depict regions where both fishing effort and the proportion of threatened species are low.

Table 1 .

 1 Species/genera selected for this study and the corresponding frequency of occurrence (percentage of trawl hauls where the taxon was caught) for the period 1994-2021. Species are sorted in descending order of occurrence. The IUCN Status column indicates the conservation status of each species: EN (Endangered), VU (Vulnerable), NT (Near Threatened), and LC (Least Concern). In this study, species are categorized as nonthreatened (0) if classified as LC or NT and threatened (1) if classified as EN or VU.

	Species	Common name	Frequency	of	IUCN	IUCN	Status	Threatened
				occurrence (%)		Status	reference
	Scyliorhinus	Small-spotted	10,7		LC	Finucci et al. (2021) 0
	canicula	catshark
	Galeus	Blackmouth		10,0		LC	Abella et al. (2016)	0
	melastomus	catshark
	Squalus	Spurdog		8,7		EN	Ellis, Soldo, et al.	1
	acanthias						(2016)
	Etmopterus	Velvet	belly	8,2		LC	Guallart et al. (2016) 0
	spinax	lanternshark
	Raja clavata	Thornback skate	7,8		NT	Ellis, Dulvy, et al.	0
							(2016)
	Chimaera	Rabbitfish		6,2		NT	Dagit & Hareide	0
	monstrosa						(2016)
	Raja asterias	Starry skate		3,6		NT	Serena, Abella, et al.	0
							(2016)
	Myliobatis	Common eagle	3,5		VU	Serena, Holtzhausen,	1
	aquila	ray					et al. (2016)
	Raja miraletus	Brown skate		3,5		LC	Dulvy et al. (2020)	0
	Mustelus spp.	Smooth-hound	3,1		VU	Farrell & Dulvy	1
							(2016)

Acknowledgments

We thank the MEDITS staff involved in the scientific sampling and analysis of biological data. This study has been conceived under the International PhD Program "Innovative Technologies and Sustainable Use of Mediterranean Sea Fishery and Biological Resources" (http://www.fishmed-phd.org/).