Conservation Status Assessment of Demersal Elasmobranchs in the Balearic Islands (Western Mediterranean) over the Last Two Decades

: More than half of the Mediterranean sharks and rays are threatened by ﬁshing exploitation. However, population assessments are limited by the scarcity of speciﬁc data on ﬁshing catches. In this study, we assessed temporal trends of the indicators developed within the European Marine Strategy Framework Directive over the last two decades in order to assess the conservation status of demersal sharks and batoids in the Balearic Islands, which represent an important fraction of the bycatch of bottom trawling in this area. On the basis of a georeferenced, ﬁshery-independent dataset of 19 species of elasmobranchs, we analyzed 20 year time series (2002–2021) of nine indicators regarding area distribution, population size, population status, and community structure. Between 30% and 50% of the elasmobranch species and functional groups showed increasing trends in distribution area and population size. This was especially true for batoids, whereas the distribution area and population size of most sharks remained stable over the study period. The remaining indicators showed stability or, in some cases, variable trends. Only in one case did we ﬁnd a negative trend sustained all along the time series (i.e., the proportion of R. radula large individuals in relation to the reference period). Overall, our results suggest that the populations of elasmobranchs from the Balearic Islands show stable or recovery trends, mainly in terms of distribution area and density. However, it remains elusive whether this community can recover to the levels of more than half a century ago, before the development of the bottom trawl ﬁshery, or whether this apparent current steady state should be interpreted as a new equilibrium.


Introduction
The vital traits (i.e., slow growth, late sexual maturation, large body size and low fecundity) of elasmobranchs (sharks and batoids) make them especially vulnerable to fishing pressure [1].In fact, fisheries are likely the main current anthropogenic threat to elasmobranch populations [2].During the 20th century, the global increase in fishing effort [3,4] led to declines and serial depletions of elasmobranch populations [5,6], which may take decades to recover [7].In the Mediterranean, more than half of the species of elasmobranchs are currently threatened [2].Recently, the General Fisheries Commission for the Mediterranean (GFCM) adopted the recommendation GFCM/42/2018/2 on fisheries management measures for the conservation of sharks and rays and the recommendation GFCM/44/2021/16 on additional mitigation measures for the conservation of elasmobranchs in the Mediterranean, including the adoption of measures to reduce their mortality by incidental catch during fishing operations, as well as supporting data collection, monitoring, and research programs on these vulnerable species [8].During the last decades, fisheries assessment and management has evolved from a species-to an ecosystem-based approach, which considers not only the fisheries target species but also other components of the ecosystem [9].This is particularly required in the complex multi-gear and multispecies fisheries of the Mediterranean, where trawling is the most important fishery in terms of fleet and catches [10].Bottom trawlers operate over a wide bathymetric range along the continental shelf and slope, causing a large level of bycatch and discards.Demersal elasmobranchs can represent an important fraction of this bycatch [11], but a minor portion of landings [12].
However, the assessment of elasmobranch populations is limited by the lack of reliable specific data on fishing catches and landings [13]; consequently, most of these populations remain unassessed and without proper management measures.As an example, a recent study by Juan-Jordà et al. (2022) [14] provided evidence of recovery of tunas and billfishes after 2000s, when fisheries management actions were implemented for target species, whereas sharks' extinction risk continues to rise due to the lack of management.In fact, the conservation status of more than 20% of sharks and rays in the Mediterranean has worsened during the last decades [15].Despite all the limitations, the use of fisheries-independent data and population and community indicators has allowed to describe the community composition and assess the status for demersal elasmobranchs in the Balearic Islands (western Mediterranean) [12,[16][17][18].These studies showed a greater diversity and density of demersal elasmobranchs in the Balearic archipelago compared to adjacent waters of the western Mediterranean, which could be related to a relatively low intensity of bottom trawling.In the Balearic Islands, the highest development of the bottom trawl fleet took place from mid-1950s to 1980, resulting in a clear decrease in elasmobranch populations and important community changes [12].During the last decades, the community of elasmobranchs has achieved some stability in terms of diversity, biomass, and specific density, but with different behavior between continental shelf and slope populations [12].Some species such as Scyliorhinus canicula, Raja clavata, and Galeus melastomus, with a wide distribution depth range, showed signs of recovery, but negative trends have been observed for other species such as Etmopterus spinax and Dipturus oxyrinchus that are more restricted to deep waters [18].These studies have related these temporal dynamics to the historical evolution of the intensity and bathymetric distribution of the bottom trawl fishing effort.
Marine biodiversity loss is considered to be one of the most severe global environmental problems.Among many other aspects, it plays an important role in maintaining marine ecosystem functioning and in providing ecosystem services.The Marine Strategy Framework Directive (MSFD; Directive 2008/56/EU) aims to increase protection of the European marine environment and to achieve a good environmental status, enabling a sustainable use of marine goods and services.Therefore, it requires the application of an ecosystem-based approach to better link ecosystem components, anthropogenic pressures, and impacts on the marine environment.Member states should analyze and periodically update the status of their marine ecosystems according to criteria and methodological standards on good environmental status (GES) and specifications of methods for their monitoring and assessment (Commission Decision EU 2017/848), established from 11 qualitative descriptors for interpreting what GES means in practice, through describing what the environment will look like when GES has been achieved.MSFD Descriptor 1 refers to biodiversity and describes the achievement of GES as the maintenance of the biological diversity.Within this descriptor, analyses of the distribution and condition of habitats and species are performed every 6 years within each marine demarcation.In regard to elasmobranchs, these periodical analyses only include the most common species at a regional level.Hence, in the Balearic Islands, the representation of the elasmobranch community in these analyses is limited to only the three most abundant species: R. clavata, S. canicula, and G. melastomus.
In the present work, we aim to expand the MSFD analysis beyond the dominant species and to evaluate the conservation status of the elasmobranch community in the Balearic Islands over the last two decades.For this purpose, we used the GES indicators for circalittoral and bathyal fishes developed in the context of MSFD Descriptor 1.A total of nine indicators regarding the population size, status, and distribution of 19 elasmobranch species were calculated on the basis of georeferenced fishery-independent data on density, biomass, and total body length.Indicators were estimated annually at the species, functional group, and community levels, and their trends were evaluated over the last two decades (2002-2021).

Study Area
The Balearic Islands are situated in the western Mediterranean, separated from the Iberian Peninsula by a distance up to 95 nautical miles, with depths between 800 and 1800 m.The continental shelf is generally narrow (3 km width), with some exceptions such as southern Mallorca, where widths can reach up to 35 km [19].The slope on the western and southern sides is gentle (6 • average inclination), while the northern and eastern sides have an abrupt slope (16 • average inclination), with a clear shelf-break.The sediments of the shelf are mainly biogenic sands and gravels, and sandy-muddy and detrital bottoms are present at the shelf-slope break, whereas muddy sediments of biogenic origin dominate the deeper areas [19,20].
The archipelago delimits the Balearic sub-basin in the north from the Algerian subbasin in the south.These sub-basins are characterized by different oceanographic conditions [21] and are connected by a series of channels with depths between 100 and 800 m, which play an important role in the regional circulation, as passages for the exchange of water masses between them [22].The Balearic sub-basin is more influenced by atmospheric forcing and Mediterranean waters, which are colder and more saline, whereas the Algerian sub-basin is affected basically by density gradients and receives warmer and less saline Atlantic waters, entering the Mediterranean through the Strait of Gibraltar [23].
Within the general oligotrophic environment of the Mediterranean, the waters around the Balearic Islands show more pronounced oligotrophy than the adjacent waters off the Iberian coast and the Gulf of Lion [24,25], due to the scarcity of rain and the karstic nature of the islands.This pronounced oligotrophy explains the high transparency of the waters in the area that enables algal populations to grow until 90-100 m depth [26,27].Red algae beds, with two main communities being rhodoliths and Peyssonnelia beds, dominate the continental shelf landscape down to 85 m depth [28,29].However, frontal mesoscale events between Mediterranean and Atlantic waters [23] and input of cold northern water into the channels [30] can act as external fertilization mechanisms that enhance productivity off the Balearic Islands.
Regarding fisheries, the Balearic Islands can be considered an individualized area for assessment and management purposes in the western Mediterranean [31], where the number of bottom trawl fishing boats has remained historically very low, compared to other areas of the Mediterranean coast off the Iberian Peninsula.Due to this lower fishing exploitation, the demersal resources and ecosystems in the Balearic Islands are in a healthier state than in adjacent areas off Iberian Peninsula, which is reflected in the population structure of the main commercial species (populations from the Balearic Islands have larger modal sizes and lower percentages of small-sized individuals), as well as in the higher abundance and diversity of elasmobranch assemblages [31].

Data Source
The data used come from the bottom trawl surveys BALAR (2002)(2003)(2004)(2005)(2006) and MEDITS (2007-2021) developed around Mallorca and Menorca in the Balearic Islands (western Mediterranean; Figure 1).These surveys are performed annually during late spring and early summer, using the experimental bottom trawl GOC-73, equipped with 10-20 mm mesh size in the cod-end, and applying a stratified sampling scheme with four depth strata (B: 51-100 m, C: 101-200 m, D: 201-500 m, E: 501-800 m) following the MEDITS protocol [32,33].Hauls were conducted during daylight hours at a towing speed between For each haul, position, depth, and the arrival and departure of the net to the bottom, as well as its horizontal and vertical openings (ranging 16-22 and 2.5-3 m, respectively), were monitored in real time using SCANMAR or MARPORT systems.Once the catch was on board, it was sorted by species, counted, and weighted, and the individual total length was measured to the lowest half cm.Abundance and biomass were standardized to 1 km 2 , using the horizontal opening of the net and the distance covered by the net on the seafloor during each haul, which was measured using a GPS.
The data used come from the bottom trawl surveys BALAR (2002)(2003)(2004)(2005)(2006) and MEDITS (2007-2021) developed around Mallorca and Menorca in the Balearic Islands (western Mediterranean; Figure 1).These surveys are performed annually during late spring and early summer, using the experimental bottom trawl GOC-73, equipped with 10-20 mm mesh size in the cod-end, and applying a stratified sampling scheme with four depth strata (B: 51-100 m, C: 101-200 m, D: 201-500 m, E: 501-800 m) following the MEDITS protocol [32,33].Hauls were conducted during daylight hours at a towing speed between 2.5 and 3 knots, with an effective duration of 20 min at depths shallower than 101 m, 30 min at depths between 101 and 200 m, and 60 min in deeper areas.
For each haul, position, depth, and the arrival and departure of the net to the bottom, as well as its horizontal and vertical openings (ranging 16-22 and 2.5-3 m, respectively), were monitored in real time using SCANMAR or MARPORT systems.Once the catch was on board, it was sorted by species, counted, and weighted, and the individual total length was measured to the lowest half cm.Abundance and biomass were standardized to 1 km 2 , using the horizontal opening of the net and the distance covered by the net on the seafloor during each haul, which was measured using a GPS.

Species Selection and Functional Groups
A total of 943 hauls from 2002 to 2021, distributed among 50 sampling stations each survey (Figure 1), were analyzed.Only elasmobranch species appearing in more than 33% of the years (at least 7 years) were considered (Table 1).Species were grouped based on

Species Selection and Functional Groups
A total of 943 hauls from 2002 to 2021, distributed among 50 sampling stations each survey (Figure 1), were analyzed.Only elasmobranch species appearing in more than 33% of the years (at least 7 years) were considered (Table 1).Species were grouped based on functional similarities, regarding their biological traits.The traits considered were type of reproduction, body shape, maximum body size, and mean body weight.Each trait was divided into four categories, and each species was assigned to one category according to Farriols et al. (2017) [34].Then, a functional resemblance matrix among species was calculated using the simple matching coefficient: where a is the number of common categories to species i and j, b is the number possessed by i and not j, c is the number possessed by j and not I, and d is the number possessed by neither.The functional resemblance matrix was subjected to cluster analysis, and a hierarchical cluster was produced using the average group linkage with PRIMER 7 [35].

Estimation of Indicators
The MSFD defines the GES for demersal fish, according to a set of criteria at the population and community levels.More related information can be found in the introduction and Article 8 of the MSFD, Directive 2008/56/EU and in the Commission Decision EU 2017/848.Here, up to nine indicators were estimated annually, six of them regarding populations and three regarding community, with the latter considered by the MSFD as ecosystem structure indicators.The meaning of each of these indicators and their calculation is detailed below.GES population indicators responding to the criteria of distribution, size, and status of populations were estimated at the species and functional group levels.Community indicators responding to the criteria of ecosystem structure were estimated: (i) for the whole community (all species included); (ii) for each functional group; (iii) for the whole community, but excluding the dominant species S. canicula, G. melastomus, and R. clavata, because their high abundance and biomass may have masked trends for the rest of the community (referred to as the reduced community hereafter).

Area of Distribution
The extent of the study area was delimited manually along the outer side of the 800 m isobath around Mallorca and Menorca and divided into 1 × 1 km cells, applying a grid established by the Spanish Ministry for the Ecological Transition and the Demographic Challenge, the authority for the implementation of the MSFD in Spain [36].
The species bathymetric range was defined by the minimum and maximum depth at which the species was caught over the whole time series.
The area of distribution indicator estimates the standardized percentage of grids where a species occurred within its bathymetric range, as follows: , where C +i is the number of grids where the species was present in year i, C t is the total number of grids within the species bathymetric range in the whole area, C mi is the number of sampled grids within the species bathymetric range in year i, and C max is the maximum number of grids within the species bathymetric range sampled any year.

Population Size Indicators Density and Biomass
Population size by year was assessed by estimating the average stratified density and biomass of each species, multiplying the average standardized density (number of individuals per km 2 ) and biomass (weight in g per km 2 ) at each stratum by the stratum surface, and the result was divided by the total surface, considering all strata.Subsequently, the annual standard score (Z-score) was calculated as follows: where Z is the Z-score value, x is the value being evaluated, µ is the mean, and σ is the standard deviation.The Z-score value for the last year indicates the current status of the species (Z last ) over the time series.

95th Percentile of Size Distribution (P95)
It is assumed that the health of a population increases as its size distribution consists of larger fish.Hence, the 95% percentile of the population length frequency distribution (P95) is considered a good state indicator for the implementation of MSFD.It is expected to be sensitive to fishing and other human impacts [37][38][39].A decrease in P95 length value generally indicates increased fishing pressure [39].
The standardized density (individuals/km 2 ) by length class was used to estimate this indicator at the species level.Species-specific size values were ordered (from smallest to largest), and the annual P95 was obtained from the length value corresponding to the ordinal (n) calculated as follows: where N i is the number of individuals in year i.

Percentage of Mature Specimens (%Mature)
For each species, the annual abundance of mature specimens was calculated as the sum of standardized density (individuals/km 2 ) of all length classes larger than the species-specific length at first maturity (L 50 ), obtained from the available literature (Supplementary Table S1).The %Mature was calculated at a species level as the ratio between the abundance of mature specimens and the total standardized density of the species.When data on sex-specific L 50 were available, this indicator was also evaluated separately for each sex.

Ecosystem Structure Indicators 2.8.1. Mean Maximum Length of Fish (MMLF)
At the community level, the proportion of large species is also a sign of good health status.The MMLF is a good indicator to reflect the evolution of the mean maximum length of the fish communities, taking into consideration intraspecific variations of size and its sensitivity to fishing pressure [39].The MMLF was calculated as follows: where L max,j is the maximum total length observed for species j, N j is the number of specimens of species j, and N is the total number of specimens.

Large Fish Indicator (LFI)
This indicator estimates the proportion of fish (relative density) with a length exceeding the length criteria to be considered a large fish.The threshold size (LFTS) that defines what could be considered a large fish within the community of elasmobranchs was determined on the basis of the proportion of individuals larger than a specific size.Different threshold sizes were tested every 5 cm from 30 to 115 cm, and the percentage of individuals larger than each threshold size was calculated.The selected large fish threshold size was that which was surpassed only by ≤2.5% of the individuals (Table 2).

Group LFTS (cm) %
Whole This indicator measures the relative abundance of individuals larger than 0.5 times the species-specific asymptotic length (L∞) in the community in relation to an initial reference period and only considering large species (i.e., those species with total maximum length larger than the large fish threshold at community level, which was here estimated at 60 cm).The asymptotic length should be calculated from the growth function by species; however, since this information was lacking in most cases, the maximum length recorded in the entire series (L max ) was used instead (L max can be consulted in Table 1).For each species, the annual abundance (stratified standardized density) of individuals larger than 0.5 times the species-specific L max was calculated.Subsequently, for each species, the mean abundance of individuals larger than 0.5 × L max in the first 3 years of the time series was computed.For the successive years, the proportion of large individuals in relation to the reference period was quantified (from here on, POL/Ref) at the species level as the ratio between the abundance of individuals larger than 0.5 × L max and mean abundance in the reference period.In cases where a species did not occur in the first 3 years, the reference period was extended until the year of first occurrence.Lastly, for any given year, the CSFb was calculated at the community level as the geometric mean of the relative abundance of the species included in the community.It should be noted that the geometric mean of a dataset with at least one zero will always be zero.Thus, given that the geometric mean does not capture any information about the non-zero values, only those species with individuals larger than 0.5 × L max occurring in all the surveys were retained for the calculation of CSFb.Therefore, the CSFb was only calculated at the whole community level considering only the species that met the aforementioned criteria: R. clavata, Raja miraletus, R. polystigma, L. naevus, Dipturus oxyrinchus, S. canicula, and G. melastomus.For all these species the reference period was the first 3 years of the time series (2002)(2003)(2004)).

Trends Evaluation
The procedure applied within the MSFD is by the means of simple linear regressions.However, taken into account that trend changes can occur throughout time series, we explored the existence of changes in the direction of the regression line, by comparing simple linear regressions with segmented models, also known as piecewise regression [40,41].In this regression analysis method, the independent variable is partitioned into intervals and a separate line segment is fit to each interval.The knot connecting two linear segments is called the breakpoint.For this purpose, the segmented R package [40,42] was used.First of all, the number of breakpoints that best fit the data was selected using the selgmented function [43] and setting at three the maximum number of breakpoints to be tested (K max ).This function performs sequential comparisons of the simple linear model using segmented models with n breakpoints (from 0 to K max ) and determines the most appropriate number of breakpoints, according to the Bayesian information criterion (BIC).When the lowest BIC was achieved with zero breakpoints, a simple linear regression model was fit to the data.Otherwise, the segmented function was used to fit a segmented model with the selected number of breakpoints.The slope, the 95% confidence interval, and the p-value (only available for simple linear regression) of the linear regression or of each linear segment were obtained.According to the sign of the slope, the range of the 95% confidence interval, and the p-value of the regression, the historical trend of each indicator was classified into three categories: increasing, decreasing, and uncertain (Table 3).For segmented models, a trend category was assigned to each interval of years.The achievement of GES was assessed on the basis of temporal trends of the indicators from 2002 to 2021.Species showing increasing or uncertain trends were considered to be in GES, while GES was not achieved by species showing decreasing trends.For segmented models, the achievement of GES was based on the last segment trend.For population size indicators (abundance and biomass), the achievement of GES was assessed on the basis of population size temporal trends and comparing the population Z last value to the mean of the annual Z-score values over the timeseries.Following the MSFD definitions of GES, the mean Z-score was assumed to be 0, and the standard deviation was assumed to be 1.Assuming that all elasmobranchs are considered K-strategists, populations were considered to be in GES when (a) the population size indicators showed a decreasing trend and Z last ≥ 0.5, or (b) the population size indicators showed an increasing or stable trend and Z last > −0.5.

Results
The elasmobranch community was clearly dominated by two sharks (S. canicula and G. melastomus) and one batoid (R. clavata) (Table 1), with different relevance depending on the bathymetric strata (Supplementary Table S2).Between 50 and 100 m depth (stratum B), S. canicula accounted for more than 85% and 54% of the density (n/km 2 ) and biomass (kg/km 2 ), respectively.In this stratum, the most abundant batoids were R. miraletus (3.6% and 2.7% of density and biomass, respectively) and R. clavate (2.6% and 16% of density and biomass, respectively).Between 101 and 200 m depth (stratum C), the community was dominated by S. canicula (accounting for 81% of density and 31% of biomass) and R. clavata (accounting for 9.4% and 40% of density and biomass, respectively).Between 201 and 500 m depth (stratum D), the most abundant species was G. melastomus, followed by S. canicula, representing 70% and 26% of density, respectively.In terms of biomass, the community of this stratum was dominated by R. clavata, G. melastmous, and S. canicula, accounting for 35%, 24%, and 24%, respectively.G. melastomus clearly dominated between 501 and 800 m depth (stratum E), representing 92% of density and 89% of biomass, followed by E. spinax (accounting for 8% of density and 4.5% of biomass) and Centrophorus uyato, with values of 0.15% and 2.7%, respectively.
The cluster analysis showed that species grouped in six functional groups (Table 1).For most species and functional groups, simple linear regressions were the models that best fitted time series of indicators (Figures 2 and 4), although there are some species and functional groups whose indicators trends were better described by segmented models.For more detail, see Supplementary Tables S3-S11.
At the functional group level, the area of distribution, stratified density, and biomass of batoids (B1, B2 and B3) increased over the study period.In contrast, the distribution and population size indicators of sharks showed apparent stability, with a few exceptions including the area of distribution of group S1, which decreased from 2002 to 2004, and the biomass of group S2, which decreased from 2010 to 2017 and increased from 2017 to 2021 (Figure 2).With regard to the population status and ecosystem structure indicators, most batoids and shark groups showed stability (Figure 2).Some groups showed different trends along the time series.This was the case of S1, which showed a decrease in P95 from 2004 and an increase in MMLF during the last years (2019-2021) and B2, along with some variability on %Mat trend during the first years (2002)(2003)(2004)(2005).At the community level, ecosystem structure indicators showed general stability (Figure 3).Considering the community as a whole, MMLF showed an increasing trend during the last part of the time series (2018-2021), a similar behavior to the MMLF of group S1.In fact, this trend was not detected when the reduced community was considered by excluding the most abundant species S. canicula, G. melastomus, and R. clavata.
At the population level, different trends were detected for group S1 before and after 2004 (Figure 4).During the first years (2002)(2003)(2004), the abundance of G. melastomus decreased, while the remaining indicators showed stability.After 2004, the distribution area of G. melastomus was stabilized, while that for S. canicula abundance increased, and P95 changed from stability to a decreasing trend.The Pol/Ref trend of G. melastomus increased during the last years of the time series (2017-2021).In group S2, although Mustelus mustelus showed an increasing trend in terms of distribution and abundance, most of the rest of indicators showed no clear trends.In group S3, trends were uncertain for all indicators evaluated, with the exception of E. spinax, which showed a decreasing trend in abundance during the first years (2002)(2003)(2004).
Within the batoid group B1, increasing trends dominated in terms of distribution and population size, with the exception of (i) Raja brachyura, which, from 2014 onward, showed negative trends of distribution and abundance, while its biomass maintained a positive trend throughout the entire period, and (ii) Dipturus oxyrinchus, which showed decreasing abundance for most of the period, except for the last 4 years, when the slope became positive.In relation to the temporal evolution of population status, Rostroraja alba presented an increasing trend in %Mature for most of the time series, and R. miraletus showed a growing trend in P95 during recent years (2014-2021).The POL/Ref indicator showed clear trends for Raja radula, with a negative slope throughout all time series, and for R. polystigma, with a sustained increasing trend.In group B2, Dasyatis spp.showed an increasing area of distribution and population size throughout the entire period and an increasing trend in POL/Ref during the last period of the time series (2018-2021).Regarding Mylobatis aquila, no clear trends were detected for any of the indicators.Group B, represented only by Torpedo marmorata, showed general signs of recovery with increasing trends for area of distribution (from 2009) and population size indicators.Note that the increasing POL/Ref trends for Torpedo marmorata should be interpreted cautiously because the L max used in this work (maximum length observed in the surveys) is much lower than that reported in the literature [44,45]; thus, it may not be representative of the whole population and not a good proxy of L∞.At the population level, different trends were detected for group S1 before and after 2004 (Figure 4).During the first years (2002)(2003)(2004), the abundance of G. melastomus decreased, while the remaining indicators showed stability.After 2004, the distribution area of G. melastomus was stabilized, while that for S. canicula abundance increased, and P95 changed from stability to a decreasing trend.The Pol/Ref trend of G. melastomus increased during the last years of the time series (2017-2021).In group S2, although Mustelus mustelus showed an increasing trend in terms of distribution and abundance, most of the rest of indicators showed no clear trends.In group S3, trends were uncertain for all indicators evaluated, with the exception of E. spinax, which showed a decreasing trend in abundance during the first years (2002-2004).Note that the CSFb was estimated for the whole community but considering only those species that surpassed the large fish threshold and occurred in all years (see Sections 2 and 2.8.3).
Overall, 30-50% of the species and functional groups analyzed showed an increasing area of distribution and population size (Figure 5).At the species level (Figure 5a), abundance was the indicator that showed more variability throughout the time series (35% of the species), while, at the group level (Figure 5b), population status indicators showed higher variability.At the species level (Figure 5a), GES was achieved by 79-100% of the species, depending on the indicator.Specifically, GES was not achieved for population size indicators of G. melastomus (abundance), Dalatias licha (abundance and biomass), C. uyato (abundance and biomass), and R. brachyura (abundance), nor for specific population status indicators such as P95 of S. canicula and CSFb of R. radula (Figure 4).GES was achieved for all indicators at the functional group (Figure 5b) and community levels (Figure 3), except for P95 of group S1 (Figure 2).
Over 40% of the species included in this work are recorded in the International Union for Conservation of Nature (IUCN) Red List of Threatened Species as vulnerable, endangered, or critically endangered in the Mediterranean (Figure 6a).Most of these species achieved the GES for the analyzed indicators (Figure 6b).However, the GES for population size was not achieved by 25% of the threatened species, which showed variable trends with negative slope on the last segment, and one species (R. radula) did not meet GES for CSFb.Over 40% of the species included in this work are recorded in the International Union for Conservation of Nature (IUCN) Red List of Threatened Species as vulnerable, endangered, or critically endangered in the Mediterranean (Figure 6a).Most of these species achieved the GES for the analyzed indicators (Figure 6b).However, the GES for  According to the temporal distribution of breakpoints of temporal trends, there are some years in which several trend changes of indicators co-occurred for different species and functional groups (Figure 7).For instance, between 2004 and 2006, there were changes in abundance trend for several species, coinciding with a change in the trend of P95.Specifically, in 2004, three species of sharks (G.melastomus, S. canicula, and E. spinax) and one species of batoid (D. oxyrinchus) improved their abundance trend, although, in the same year, one species suffered a negative trend change for P95.Between the years 2011 and 2012, changes in distribution, population size, P95, and %Mat occurred for different species.In 2014 and 2017-2018, some trend changes also co-occurred.All these periods of change were also detected at the functional group level.According to the temporal distribution of breakpoints of temporal trends, there are some years in which several trend changes of indicators co-occurred for different species and functional groups (Figure 7).For instance, between 2004 and 2006, there were changes in abundance trend for several species, coinciding with a change in the trend of P95.Specifically, in 2004, three species of sharks (G.melastomus, S. canicula, and E. spinax) and one species of batoid (D. oxyrinchus) improved their abundance trend, although, in the same year, one species suffered a negative trend change for P95.Between the years 2011 and 2012, changes in distribution, population size, P95, and %Mat occurred for different species.In 2014 and 2017-2018, some trend changes also co-occurred.All these periods of change were also detected at the functional group level.According to the temporal distribution of breakpoints of temporal trends, there are some years in which several trend changes of indicators co-occurred for different species and functional groups (Figure 7).For instance, between 2004 and 2006, there were changes in abundance trend for several species, coinciding with a change in the trend of P95.Specifically, in 2004, three species of sharks (G.melastomus, S. canicula, and E. spinax) and one species of batoid (D. oxyrinchus) improved their abundance trend, although, in the same year, one species suffered a negative trend change for P95.Between the years 2011 and 2012, changes in distribution, population size, P95, and %Mat occurred for different species.In 2014 and 2017-2018, some trend changes also co-occurred.All these periods of change were also detected at the functional group level.95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation to the reference period.CSFb: conservation status of fish.The size of the circles shows the number of trend changes (breakpoints) at each year.Each color refers to one indicator, detailed in the y axis.

Discussion
In a highly multispecies fishery such as Mediterranean bottom trawling [10], the conservation status of the bycatch species should be monitored and included in the fishery assessment and management.This is especially relevant within the current context of the Ecosystem Approach to Fisheries [9] and even more in the case of vulnerable species, such as sharks and batoids, which represent an important fraction of the bycatch of this fishery [11,46,47].
Despite the high vulnerability and worsening threat status of elasmobranchs [14], this is a poorly studied group.The lack of reliable catch time series makes it not feasible to assess the status of their populations through stock assessment models, as conducted for target species of teleosts, decapod crustaceans, and cephalopods.However, fishery-independent scientific surveys provide valuable standardized data of nontarget species such as elasmobranchs.Here, we applied the indicators developed in the frame of the MSFD, to assess the conservation status of the demersal elasmobranch community in the Balearic Islands beyond its dominant species using data from the MEDITS survey, a program of scientific surveys developed with bottom trawl along the whole northern Mediterranean and included in the European framework for the collection, management, and use of data in the fisheries sector to support scientific advice to the Common Fisheries Policy [32].
In general, the elasmobranch community and populations assessed here showed nonsignificant trends for most of the analyzed indicators, and even showed some signs of recovery for several species, mostly batoids and mainly in relation to population distribution and size.Sharks did not show clear dominant trends, except for some specific cases, such as G. melastomus, S. canicula, and E. spinax, with negative slopes at the beginning of the time series, which, after 2004-2005, stabilized or even increased trends for the first two shark species.Regarding population status indicators, trends were not significant except for some particular cases.
Significant trends were mainly detected in the species with the highest frequency of occurrence (e.g., sharks from group S1 and most batoids from group B1), with some exceptions of less common species (e.g., R. alba), especially with regard to the indicators of distribution and population size.In the initial National Assessment Report for the implementation of the MSFD in Spain, only three elasmobranchs, the most common S. canicula, G. melastomus, and R. clavata, were included in Descriptor 1 for the Eastern and Balearic Demarcation, which covers the northeastern Iberian Peninsula and the Balearic Islands [48].For each indicator, population temporal trends were described and classified into three categories: increasing, decreasing, and stable.The first two categories referred to significant trends with positive and negative slopes, respectively, while the latter category was assigned to species showing nonsignificant trends.In the present study, we expanded the MSFD evaluation beyond the most common species, including those with lower frequency of occurrence for which data are not available for all years (e.g., C. uyato and D. licha).For this reason, and to be more conservative regarding the potential scientific advice for the management of these vulnerable species, on the basis of our results, we adapted trend categories, classifying the nonsignificant trends as uncertain instead of stable.In the case of the most common species with nonsignificant slopes, herein classified as uncertain, trends could be considered stable, whereas, for less frequent species, the trend could not be described with certainty.
Elasmobranchs, due to their K-selected life-history strategies and their high trophic position [49], are species particularly vulnerable to fisheries and can be considered indicators of fishing pressure.Although it was already known that Balearic Islands have higher diversity and density of the demersal elasmobranch community than the adjacent waters of the Iberian Peninsula, probably due to a lower intensity in fishing exploitation together with some biogeographic factors [16,17], the apparent general stability or recovery trends described here contrast with the global and Mediterranean trends [2,14,15].This may be explained by the temporal evolution of the fishing effort in this area, especially regarding the bottom trawl fleet, which represents the most important demersal fishery around the Balearic Archipelago [31].This fleet uses low selective fishing gears, with high bycatch and discard rates, thus having a high impact on the benthic and nektobenthic species, including elasmobranchs [11,50,51].
In this area, despite the sustained decrease in the number of bottom trawlers since the mid-1970s, the fishing capacity remained quite constant, and even higher, due to the continuous increase in the engine power of the vessels [52].In fact, their landings did not show any decrease, and they remained fairly stable.However, during the last decade, the engine power of the vessels did not increase, mainly due to the continuous increase in oil prices (up to 45%), which represents a high percentage of running costs of the trawl fleet, but the number of vessels continued to decline [52].Indeed, the analysis of the Vessel Monitoring by Satellite System (VMS), since its implementation in this fleet in 2006, showed a constant decrease of the number of trawl fishing days around the Balearic Islands (Figure 8).This reduction, which represents up to 60%, is especially pronounced in the intermediate slope fishing grounds (Figure 8, strata E), between 501 and 800 m depth.This could be explained by the greater distance between these fishing grounds and the ports and some episodes of sudden and strong reductions in the catches of red shrimp (Aristeus antennatus), the target species of bottom trawling at this depth, mainly due to environmental factors influencing the availability of this species to fishing exploitation [53,54], which have occurred in the study area (personal observations).Moreover, because the economic crisis during the recent years, consumption habits could have changed, reducing the demand for deep water decapod crustaceans of high economic value and increasing the demand for cheaper species, fish and cephalopods, which are caught on the continental shelf [55,56].The continental shelf and upper slope (strata B-D, between 51 and 500 m depth) showed a fairly stable effort trend, but also a slightly decreasing one after 2013-2014.Therefore, it would be reasonable to attribute the observed trends here described in the elasmobranch community to the decrease in trawl fishing effort.Our results agree with previous studies developed in the same area, but from different approaches [12,18], which also suggest similar optimistic recent trends for the elasmobranch community from the continental shelf, but not for deep water species.These authors attributed the increase in some slope-dwelling species during the last decades to Our results agree with previous studies developed in the same area, but from different approaches [12,18], which also suggest similar optimistic recent trends for the elasmobranch community from the continental shelf, but not for deep water species.These authors attributed the increase in some slope-dwelling species during the last decades to (i) the change in the bathymetric distribution of the fishing effort due to the bottom trawl fleet displacement toward greater depths, targeting the red shrimp, which has resulted in an effort reduction on the continental shelf bottoms [57,58], and (ii) different degrees of resilience of species [1].
Environmental factors could also be playing a role in determining population dynamics and trends.It is expected that the multiage population structure of elasmobranchs could buffer to some extent the effects of environmental fluctuations [59].However, fishery exploitation can lead to truncated age structures, increasing population sensitivity to environmental factors [60].Climate change has been linked to the distributional shifts [61,62] and bathymetric range reductions [63] of elasmobranchs.
Despite the increasing abundance trend of S. canicula, our results showed a decrease in the 95th percentile of its size frequency distribution (P95) after 2004.A decrease in the value of size indicators is usually a sign of an increase in fishing pressure [39]; however, as previously mentioned, the bottom trawl fishing effort, the fishery that exploits this species almost exclusively in the Balearic Islands, has decreased since 2011.Interestingly, a decreasing trend in the length of first maturity for S. canicula was detected in the area, suggesting an evolutionary response to overfishing [18].Fishing mortality can lead to changes in population dynamics, through compensatory changes or evolutionary response to overfishing [1].Increased survival of juveniles, rather than increased fecundity, also provides greater resilience to fishing pressure [64]; thus, changes in P95 could be either a direct consequence of the fishery removal of larger specimens or an evolutionary response to it.
According to the temporal distribution of breakpoints, some periods of generalized changes in elasmobranch populations trends were identified.In 2004-2005, several trend changes co-occurred.In agreement, Ramirez-Amaro et al. (2020) [18] also described relevant changes in the diversity indices until 2004, when trends became stable.Between 2011 and 2012, several trend changes also took place, in this case coinciding with the decrease in the trawl fishery total effort.
The populations of threatened species showed uncertainty in general.Increasing trends were described for 30% of these species, e.g., M. mustelus or R. alba.Some species showed decreasing trends (e.g., R. radula for POL/Ref) or did not achieve GES for the population size indicators, e.g., the critically endangered C. uyato or the vulnerable D. licha.The GFCM recommendations require high protection from fishing activities for elasmobranch species listed in Annex II of the Barcelona Convention (GFCM/42/2018/2), as in the case of R. alba, including species-specific actions for some species (GFCM/44/2021/16).Given that, a bigger percentage of increasing trends would be desirable, and efforts should be made to improve the conservation status of threatened populations.
The Mediterranean has suffered intense fishing for centuries [65].The community composition described here matches the descriptions from other recent studies [18], but it shows major changes compared to the elasmobranch community before the development of industrial fisheries in the 1960s.Species such as S. canicula and G. melastomus, which currently predominate in the elasmobranch assemblages and represent a significant fraction of their current catches, may have had lower occurrence and abundance.Bottom trawling can change benthic communities and modify the abundance of the preferred preys of these species, as well as increase the availability of fishing discards on which they have also been observed to feed in the northeast Atlantic [66][67][68].This could also happen in the Mediterranean, where a very generalist diet for S. canicula [69] and a change in the feeding habits of G. melastomus have been suggested, due to the additional energy provided by fishing discards [70,71].In contrast, previously fairly common species, such as Scyliorhinus stellaris, Galeorhinus galeus, and Oxynotus centrina have practically disappeared from current fishing catches, while other species such as Torpedo torpedo, Squalus acanthias, Rhinobathos spp., and Squatina spp.can even be considered locally extinct or near extinction [72].Thus, the current apparent steady state could be interpreted as a new equilibrium.In this sense, it should be noted that shifting baselines may lead to failure to perceive true levels of decline [73].Historical records and studies could provide more realistic perspective for the interpretation of current population status.In the northwestern Mediterranean, historical data have shown a clear decline in demersal elasmobranch species populations since the 1960s [74].Given that, trends reported in studies covering only the last decades should be interpreted carefully.For instance, according to Aldebert, M. mustelus suffered an important decrease at the beginning of the intense industrial fishing period, due to its economic value; thus, the increasing trend in distribution and population size reported in the present work should not be interpreted as evidence of good conservation status, but a sign of possible recovery if the trend is maintained into the future.
The use of a set of indicators developed within the framework of the MDSF allowed us to obtain an evaluation of the status of conservation of the elasmobranch community, at both a general and a species level, using agreed parameters to define the GES.Given that elasmobranchs are K-strategists, the GES definitions for vulnerable species have been applied, despite not all elasmobranchs being equally vulnerable.Future evaluations could take into consideration the lower potential of recovery from exploitation of large and/or late-maturing species [15] such as D. oxyrinchus, R. alba, and E. spinax [75][76][77].In the MSFD frame, the GES is defined on the basis of indicator trends and, when applicable, in relation to threshold values.In this study, as well as in the initial evaluation for the implementation of the MSFD [48], reference values were used for the indicators of population size and fish conservation status.Defining GES in this way gives interesting results since, for example, GES was not achieved by G. melastomus, D. licha, and C. uyato in terms of population size indicators, despite showing increasing or uncertain trends.Further work is required to define threshold values for all indicators.

Conclusions
Although elasmobranch populations in the Mediterranean have generally decreased over the last decades, the demersal sharks and rays of the Balearic Islands present certain stability and even signs of recovery, probably due to the reducing bottom trawl fishing intensity in the area.This was observed with regard to the distribution and size of batoid populations.For the indicators of the state of the population and their conservation, no clear trends were found, except for specific, particularly relevant cases.This is the case for species cataloged as vulnerable or endangered, such as M. mustelus or R. alba, respectively, which show clear increasing trends for the indicators of abundance and percentage of mature individuals, respectively.Altogether, these results allow for some optimism about the sustainability of elasmobranch populations around the Balearic Islands and highlight the importance of monitoring and assessing bycatch species in order to generate the required scientific knowledge to develop and apply the Ecosystem Approach to Fisheries.However, it remains difficult to predict whether the demersal elasmobranch community of the Balearic Islands will recover to the levels of more than half a century ago or whether this apparent current steady state should be interpreted as a new equilibrium.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/fishes8050230/s1,Table S1: Summary of information related to length at first maturity (L50) for demersal elasmobranchs distributed in the Balearic Islands; Table S2: Total standardized density and standardized biomass for each stratum; Table S3: Density temporal trend fitted through linear or segmented regressions (slope, p-value, 95% confidence interval, and adjusted R 2 ); Table S4: Density temporal trend fitted through linear or segmented regressions (slope, p-value, 95% confidence interval, and adjusted R 2 ); Table S5: Area of distribution temporal trend fitted through linear or segmented regressions (slope, p-value, 95% confidence interval, and adjusted R 2 ); Table S6: The 95th percentile temporal trend fitted through linear or segmented regressions (slope, p-value, 95% confidence interval, and adjusted R 2 ); Table S7: Temporal trend of percentage of mature specimens fitted through linear or segmented regressions (slope, p-value, 95% confidence interval,

Fishes
3 knots, with an effective duration of 20 min at depths shallower than 101 m, 30 min at depths between 101 and 200 m, and 60 min in deeper areas.

Figure 1 .
Figure 1.Location of the study area (Mallorca and Menorca, Balearic Islands, western Mediterranean), showing the BALAR/MEDITS stations (blue dots) from 2002 to 2021.The isobaths of 50, 100, 200, 500, and 800 m (lines in the periphery of islands), corresponding to MEDITS strata depth ranges, are shown.

Figure 1 .
Figure 1.Location of the study area (Mallorca and Menorca, Balearic Islands, western Mediterranean), showing the BALAR/MEDITS stations (blue dots) from 2002 to 2021.The isobaths of 50, 100, 200, 500, and 800 m (lines in the periphery of islands), corresponding to MEDITS strata depth ranges, are shown.

Figure 2 .
Figure 2.Trends (increasing: upward green arrow, decreasing: downward red arrow, uncertain: flat yellow arrow) and good environmental status achievement (GES: Y, achieved; N, not achieved) at the group (S1: sharks 1, S2: sharks 2, S3: sharks 3, B1: batoids 1, B2: batoids 2, B•: batoids 3) and community levels.Whole community: all species considered.Reduced community: excluding the dominant species Scyliorhinus canicula, Galeus melastomus, and Raja clavata.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the speciesspecific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.CSFb: conservation status of fish.For population size indicators (density and biomass), the achievement of GES was assessed on the basis of population size temporal trends and comparing the population Zlast value to the mean of the annual Z-score values over the time series as follows: populations were considered to be in GES when (a) the population size indicator showed a decreasing trend and Zlast ≥ 0.5, or (b) the population size indicator showed an increasing or stable trend and Zlast > −0.5.For the remaining indicators (distribution and population status), only the temporal trends were taken into consideration (groups showing increasing or uncertain trends were considered to be in GES).

Figure 2 .
Figure 2.Trends (increasing: upward green arrow, decreasing: downward red arrow, uncertain: flat yellow arrow) and good environmental status achievement (GES: Y, achieved; N, not achieved) at the group (S1: sharks 1, S2: sharks 2, S3: sharks 3, B1: batoids 1, B2: batoids 2, B3: batoids 3) and community levels.Whole community: all species considered.Reduced community: excluding the dominant species Scyliorhinus canicula, Galeus melastomus, and Raja clavata.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.CSFb: conservation status of fish.For population size indicators (density and biomass), the achievement of GES was assessed on the basis of population size temporal trends and comparing the population Z last value to the mean of the annual Z-score values over the time series as follows: populations were considered to be in GES when (a) the population size indicator showed a decreasing trend and Z last ≥ 0.5, or (b) the population size indicator showed an increasing or stable trend and Z last > −0.5.For the remaining indicators (distribution and population status), only the temporal trends were taken into consideration (groups showing increasing or uncertain trends were considered to be in GES).

Figure 3 .
Figure 3. Trends at community level for the indicators.The first row of plots (a-c) corresponds to the whole community; the second row (d,e) corresponds to the community, excluding the dominating species (Scyliorhinus canicula, Galeus melastomus, and Raja clavata).The first column (a,d) shows trends of the large fish indicator (LFI); the second column (b,e) shows the trends of the mean maximum length of fish (MMLF); the third column (c) shows the trend of the conservation status of fish (CSFb).Note that the CSFb was estimated for the whole community but considering only those species that surpassed the large fish threshold and occurred in all years (see Materials and Methods, Section 2.8.3).

Figure 3 .
Figure 3. Trends at community level for the indicators.The first row of plots (a-c) corresponds to the whole community; the second row (d,e) corresponds to the community, excluding the dominating species (Scyliorhinus canicula, Galeus melastomus, and Raja clavata).The first column (a,d) shows trends of the large fish indicator (LFI); the second column (b,e) shows the trends of the mean maximum length of fish (MMLF); the third column (c) shows the trend of the conservation status of fish (CSFb).Note that the CSFb was estimated for the whole community but considering only those species that surpassed the large fish threshold and occurred in all years (see Sections 2 and 2.8.3).

Figure 4 .
Figure 4. Trends (increasing: upward green arrow, decreasing: downward red arrow, uncertain: flat yellow arrow) and good environmental status achievement (GES: Y, achieved; N, not achieved) at the species level.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).POL/Ref: proportion of large individuals in relation to the reference period.For population size indicators (abundance and biomass), the achievement of GES was assessed on the basis of the population size temporal trends and comparing the population Zlast value to the mean of the annual Z-score values

Figure 4 .
Figure 4. Trends (increasing: upward green arrow, decreasing: downward red arrow, uncertain: flat yellow arrow) and good environmental status achievement (GES: Y, achieved; N, not achieved) at the species level.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature

Fishes 2023, 8 , 24 Figure 5 .
Figure 5. Number of species (a) and groups (b) by trend type for each indicator and number of species (a) and groups (b) achieving the good environmental status (GES).P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the speciesspecific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation to the reference period.Trend "variable" refers to species showing trend changes throughout the time series (segmented regressions).

Figure 5 .
Figure 5. Number of species (a) and groups (b) by trend type for each indicator and number of species (a) and groups (b) achieving the good environmental status (GES).P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the speciesspecific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation to the reference period.Trend "variable" refers to species showing trend changes throughout the time series (segmented regressions).

Figure 6 .
Figure 6.Threatened species.(a) Percentage of species for each IUCN status category; (b) histogram showing the number of threatened species according to the type of temporal trend for each indicator.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation to the reference period.Trend "variable" refers to species showing trend changes throughout the time series (segmented regressions).

Figure 6 .
Figure 6.Threatened species.(a) Percentage of species for each IUCN status category; (b) histogram showing the number of threatened species according to the type of temporal trend for each indicator.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation to the reference period.Trend "variable" refers to species showing trend changes throughout the time series (segmented regressions).

Fishes 2023, 8 , 24 Figure 6 .
Figure 6.Threatened species.(a) Percentage of species for each IUCN status category; (b) histogram showing the number of threatened species according to the type of temporal trend for each indicator.P95: 95th percentile of the size frequency distribution.%Mature: percentage of mature specimens (larger than the species-specific length at first maturity, L50).MMLF: mean maximum length of fish.LFI: large fish indicator.POL/Ref: proportion of large individuals in relation the reference period.Trend "variable" refers to species showing trend changes throughout the time series (segmented regressions).

Figure 7 .
Figure 7. Temporal distribution of changes in trend (breakpoints) for the different indicators assessed at the species level (filled circles), group level (empty circles), and community level (striped).P95:

Figure 8 .
Figure 8. Temporal trends of bottom trawl fishing effort in the Balearic Islands.Number of fishing days from VMS records in total (red and wider arrow) and for each depth stratum (B: 50-100 m, C: 101-200 m, D: 201-500 m, E: 501-800 m).

Table 3 .
Categories of temporal trends and criteria used to define trends.Slope: sign of the regression slope.95% CI range: range of the 95% confidence interval.Symbol: symbol used in the results table to describe trends.* The p-value is only available for simple linear regression.