The diversity of recent trends for chondrichthyans in the Mediterranean reflects fishing exploitation and a potential evolutionary pressure towards early maturation

Chondrichthyans are a vulnerable group that has been overexploited for almost half a century in the Mediterranean. Since in this area most chondrichthyans are rarely incorporated into international statistics, the impact of fishing on their populations is difficult to assess. Here, we evaluate temporal trends in order to understand the recent history of chondrichthyans in the western Mediterranean. Fishery-independent data were obtained from scientific surveys carried out from 1994 to 2015 in three geographical sub-areas. Our results reflect fairly stable populations in terms of diversity, with some increase in density and standardized biomass of some species dwelling on the continental shelf, and even for some species dwelling on the slope. In contrast, decreasing trends were observed in some deep-water species. This can be explained by the reduction of the trawling effort on the continental shelf over the last few decades, and the shift of the fleet towards deep waters, along with the greater resilience displayed by some species. Furthermore, a decreasing trend in maturity of Scyliorhinus canicula was detected, suggesting an evolutionary response to overfishing. These results improve scientific knowledge for developing true adaptive management in the Mediterranean that will implement measures to strengthen or initiate the recovery of chondrichthyans.


Results
A total of 34 chondrichthyan species belonging to 13 families were caught throughout the western Mediterranean over the last two decades, 27, 26 and 19 species of which were caught in GSA01, GSA05 and GSA06, respectively (Table 1).
Diversity and species composition. Time series exhibited similar trends for each diversity index (S, H' and J'), with important fluctuations during the first 15 years, followed by stability over the last five years (Fig. 2,  Supplementary Figure S1). DFA results of diversity indices are given in Supplementary Table S2.
Scyliorhinus canicula and Galeus melastomus were the most important sharks, and even chondrichthyans, in all GSAs. Scyliorhinus canicula represented >90% and >70% in terms of density and standardized biomass, respectively, in bathymetric strata B and C. In stratum D, G. melastomus corresponded to >50% in both density and standardized biomass in GSA01 and GSA05, while in GSA06, the most frequent was S. canicula (54% and 82% in density and standardized biomass, respectively). In stratum E, G. melastomus represented >85% in both density and standardized biomass (Supplementary Table S3).
The most important batoids in both density and standardized biomass showed differences not only by depth strata, but also by GSA (Supplementary Table S4). In GSA01, the most frequent species were Torpedo marmorata (>40% in both density and standardized biomass) and Raja asterias (14% and 55%, respectively) in stratum B, whilst Leucoraja naevus predominated in strata C and D (>70% in both density and standardized biomass). In GSA05, the most important species in stratum B were Raja miraletus (25% and 5% in density and standardized  85 , and mean annual landings per unit effort (blue line; with standard error bars) of elasmobranchs in the Balearic Islands 24 . Colours correspond to different stages of the catches identified according to the algorithm for catch-based stock status 86 : green is "developing", yellow is "exploited" and red is "overexploited". The black square indicates our study period. values were observed in GSA05 and GSA06, suggesting stability over the study period. Stratum E exhibited cyclic changes throughout the study period. In this stratum, negative factor loadings were observed in the case of GSA01, indicating an inverse trend compared to the common trend (Fig. 3, Supplementary Figure S2). The two most frequent species of sharks, S. canicula and G. melastomus, also showed increasing common trends in both density and standardized biomass over the last 10 years for all GSAs. Nonetheless, low factor loadings observed in density and standardized biomass for GSA05 and GSA02 respectively, suggesting stability throughout the whole period analysed (Fig. 4, Supplementary Figure S3). In the case of Etmopterus spinax, DFAs revealed two possible common trends in density, and one in standardized biomass. The trend in standardized biomass and second trend in density (Fig. 4) showed a decrease over the last 10 years for GSA05. In GSA01 and GSA06, the first trend in density exhibited stability, with fluctuations throughout the study period. DFA results for sharks are given in Supplementary Tables S5-S6. In general, batoids also showed a fairly common increasing temporal trend for the three GSAs, with a rise in density and standardized biomass over the last 10 years (Fig. 5, Supplementary Figure S4). GSA01 and GSA05 showed low factor loadings in strata C and D, respectively, exhibiting stability with fluctuations throughout the study period. Torpedo marmorata showed a decrease in both density and standardized biomass over the first five www.nature.com/scientificreports www.nature.com/scientificreports/ years, followed by an increase for the next eight years, and finally stability over the last five years (Fig. 5). DFA results for batoids are given in Supplementary Tables S7-S8. In GSA05, the most frequent batoid species, R. clavata, R. miraletus, R. polystigma, R. radula, L. naevus and D. oxyrinchus, exhibited different patterns (Fig. 6). Raja clavata showed an increase over the study period in both density and standardized biomass, while R. miraletus, R. www.nature.com/scientificreports www.nature.com/scientificreports/ polystigma, R. radula and L. naevus revealed stability and D. oxyrinchus showed a decrease in density, but stability for standardized biomass.
Galeus atlanticus and Chimaera monstrosa showed stability in both density and standardized biomass (Fig. 6). The first species was only present in GSA01, while the second one was also caught regularly at this GSA.
Size. DFA models of size structure could be analysed for three species of sharks, S. canicula, G. melastomus and E. spinax. The results of DFA models of these species are given in Supplementary Tables S9-S11. In the case of batoids, the only species whose size structure could be analysed was R. clavata in GSA05.
For S. canicula, five length categories were modelled. In the category <20 cm TL, two common trends were observed. The first one was observed in GSA01 and GSA06, showing fluctuations over the early years, followed by stability in the last five years. The second one was observed in GSA05, with stability throughout the study period.  www.nature.com/scientificreports www.nature.com/scientificreports/ The same length categories were analysed for G. melastomus. The categories <20 and 20-30 cm TL showed similar trends, with fluctuations in density across the study period, and a slight increase over the last five years. The remaining categories (30-40, 40-50 and >50 cm TL), in general showed stability with slight fluctuations over www.nature.com/scientificreports www.nature.com/scientificreports/ the study period. In the last two categories, GSA05 had a low negative factor loading, indicating a stability over the study period (Fig. 7, Supplementary Figure S5).
In E. spinax, four length categories were analysed. The category <15 cm TL showed a common trend with a sharp drop followed by a slight increase until 2007 and then stability until the end of the series. The category  www.nature.com/scientificreports www.nature.com/scientificreports/ In the case of R. clavata, five length categories were analysed. The categories 40-50 and 50-60 cm TL showed increasing trends, while the remaining categories (<30, 30-40 and >60 cm TL) displayed stability over the study period (Fig. 8).
Maturity. Length at first maturity (L 50 ) could be estimated for two sharks S. canicula and G. melastomus, in all GSAs, and one batoid, R. clavata, in GSA05. Parameters of maturity ogives by year, species, and sex are given in Supplementary Tables S12-S14.
Linear regression models showed a decrease in L 50 of S. canicula for both sexes and in all GSAs (Fig. 9). In females, this parameter decreased around 4 cm TL over 16 years in GSA01 and GSA06, and around 1.5 cm TL over 9 years in GSA05. For males this reduction was around 2 cm TL across the same periods. For G. melastomus, a decrease of L 50 was also detected for females in GSA06 and males in GSA01, with a reduction of around 2 cm TL over 16 years in both sexes. By contrast, R. clavata did not show any trend in L 50 over the study period. The size of the smallest mature female showed similar trends for the three species analysed (Fig. 9).
Since S. canicula was the only species with a significant reduction in L 50 in both sexes and in the three GSAs, its somatic factor was estimated for small-sized individuals (<35 cm TL). With the only exception of females in GSA01, no significant trends were found across the study period (Fig. 10).

Discussion
In the Mediterranean, there has been increasing international concern for severe declines in the abundance and diversity of chondrichthyans, calling for urgent assessment of fishery exploitation. Here we applied a novel approach, at community and population level, to evaluate the temporal trends of demersal chondrichthyans in three geographical sub-areas along the western Mediterranean over the last two decades, using diversity, density, standardized biomass and biological parameters.
While the diversity indices used in this study measure different aspects of diversity, all of them exhibited similar trends, with relevant changes the firsts years (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) and stability the last 15 years for the three geographical sub-areas, such as has been observed in the diversity of demersal fishes 34 and cephalopods 35 throughout the whole Mediterranean. This diversity stability could be explained by the fact that relevant changes in diversity of demersal chondrichthyan communities may have occurred long before the onset of scientific monitoring of fisheries as has been reported in demersal fish communities 36 . Considering the available scientific literature in the Mediterranean 16,17,[37][38][39][40] , there are differences in the demersal chondrichthyan community composition between early fishing explotation (1950's) and the period studied here . In early fishing exploitation, some of the most threatened chondrichthyans currently such as Galeorhinus galeus, Oxynotus centrina, Scyliorhinus stellaris, Squalus acanthias, Squatina spp. and Rhinobathos spp. 18,19 were common in descriptions of fish communities in the trawl fishing grounds off the Balearic Islands 37-40 and the Adriatic Sea 17 . However, these species are no longer detected in scientific surveys for the assessment of demersal resources and have almost disappeared from commercial landings 16 , in some cases being considered locally extinct 41 . This change in chondrichthyan community can also be reflected in the percentage of species occurrence, wherein in early fishing exploitation period Scyliorhinus canicula and Galeus melastomus represented around 40% and 10% 16,17 of species occurrence and its percentage of occurrence has recently increased to around 90 and 60%, respectively.
In contrast with the decreasing population trends described for chondrichthyan species in different areas of the Mediterranean, such as the Gulf of Lion 16 , the Ligurian Sea 13 , the Adriatic Sea 4,14 and the Aegean Sea 12 , our results of density and standardized biomass suggested generally stable trends. Signs of recovery were even observed for the assemblages and the most abundant species inhabiting the continental shelf (S. canicula, Torpedo mamorata and Raja clavata) and for the most abundant species dwelling on the slope (G. melastomus). By contrast, decreasing trends were only detected for the deep water species, Etmopterus spinax and Dipturus oxyrinchus.
Fishing pressure could be one of the main factors contributing to the observed trends. The temporal evolution of the bottom trawl fleet in the western Mediterranean has gone through three phases 8,27 : (i) from the mid 1960s to mid 1970s, the number of vessels increased by a factor of 2.5; (ii) from the mid 1970s to 2000, this number decreased gradually, although the fishing capacity of this fleet remained constant or even increased, due to its continued growth in terms of engine power; and (iii) since 2000, the number of vessels has continued to decrease, probably along with the fishing capacity of the fleet, as mean engine power has remained fairly constant, especially during recent years. In this sense, our study covers a period of time when the chondrichthyan communities had already been altered in their diversity, and the effort of bottom trawl fishery was fairly stable or, at least, not increasing as previously. However, it is also necessary to consider the technological improvement of fishing gears and vessels during recent decades, which has enabled the worldwide expansion of fisheries to new areas, such as the open ocean and deep sea bottoms 42 . This has occurred in the western Mediterranean, where trawl fishery has moved from the continental shelf to the slope to capture the decapod crustaceans red shrimp (Aristeus antennatus) and Norway lobster (Nephrops norvegicus) with higher economic value 21,43 . The no increase in fishing effort of the bottom trawl, jointly with its displacement to deep waters, must be the basis for the current stability and recovery shown by sharks and batoids on the continental shelf. This shift in trawl activity has also been suggested to explain the recent stability observed in demersal chondrichthyans off Sardinia 44 , in the central western Mediterranean. Similarly, in the north-east Atlantic the progressive reduction of fishing pressure has resulted in the recovery of certain stocks, including S. canicula [45][46][47] .
The two most abundant demersal chondrichthyans in the western Mediterranean, S. canicula and G. melastomus, showed increasing trends despite their different bathymetric distribution (continental shelf and shelf break, and upper and middle slope, respectively 31,48 ). Contrary to large chondrichthyans, these two small-medium sharks show life history traits that are less susceptible to fishing pressure (e.g. early maturity, shorter generation time, fast growth, continuous reproductive cycles 49-51 ), giving them greater resilience and capacity of recovery 52 , as has been recently documented in small coastal sharks in the north-western Atlantic 53  www.nature.com/scientificreports www.nature.com/scientificreports/ and G. melastomus are considered to be opportunistic scavengers that may modify their natural diet to benefit from the disturbed sediments and discards generated by bottom trawling 54 ; moreover, the former species has shown a high survival rate as a discarded catch of this fishery 55 . A similar situation has been observed for some batoids captured on the continental shelf (R. brachyura, R. clavata, R. montagui and L. naevus 56 ), whose populations showed stability and even recovery in the present study.
The only decreasing trends were observed in the deep-water chondrichthyans, D. oxyrinchus and E. spinax. The decrease in the latter had already been reported in the north western Mediterranean 31,57 . Deep water www.nature.com/scientificreports www.nature.com/scientificreports/ chondrichthyans, with longer turnover times than the species dwelling on the continental shelf, are considered to be especially vulnerable to fishing exploitation 58,59 . This is the case of E. spinax, with late maturation, long reproductive cycle, low fecundity, and considerable longevity 60,61 . In this sense, E. spinax is quite different from G. melastomus, both species showing distinct life strategies and opposite trends. The deeper bathymetric distribution of G. melastomus, which in the western Mediterranean is relatively abundant down to 1400 m depth [62][63][64] , could also contribute to the stability of populations of this species. In this area, bottom trawl fishery does not reach more than 800-900 m depth 65 and, moreover, the GFCM decided in 2006 to ban this fishery beyond 1000 m depth.
Regarding biological descriptors of main species, no clear trends were observed in size composition, although a decreasing trend in length at first maturity (L 50 ) and the size of the smallest mature female was detected for S. canicula in all GSAs. It is known that fishing overexploitation could alter size structure and population parameters (e.g. fecundity, maturity, growth rate) on chondrichthyan populations in response to changes in species abundance (density-dependence change) 3,66-68 . Early maturity detected in S. canicula could be an indicator that the populations of this species have been under stress, due to the high level of trawl fishing exploitation in the study area. In fact, general overfishing of the demersal population around the Balearic Islands during the 1980s included elasmobranchs 8 . This early maturation, an adaptive response that increases the likelihood of offspring reaching maturity 69 , could have promoted the present recovery of S. canicula. Two possible mechanisms could explain this adaptation: an evolutionary response to selection for smaller size at maturity and/or by reducing intra-specific competition, due to reduction in species abundance, with the consequent increase in available food 66,70 . Environment is another potential mechanism that is deemed to influence life history changes in chondrichtyans 3,27 . Favourable conditions (e.g. intensive upwelling, high productivity) produce higher quality offspring that grow at faster rates, resulting in early maturation 67,68 . However, since no significant trends were found in somatic condition, this is not the case of S. canicula populations. Thus, the most plausible hypothesis to explain the early maturation of this species is its evolutionary response to overfishing. In fact, before changing growth www.nature.com/scientificreports www.nature.com/scientificreports/ rate, fishes use their somatic energy reserves to meet the running costs of their basal metabolism (e.g., maintenance, immune defence, and cognition), digestion, or routine activities such as migration to breeding grounds, courting, competition for mates and breeding sites, copulation, and parental care 71,72 .
Detecting changes in population parameters of chondrichthyans is difficult since they are typical k strategists in comparison to teleosts 3,73 . Given that small size chondrichthyans have a higher intrinsic population growth rate than larger species 52,73 , these changes have been detected mainly in few small-medium sharks such as Rhizoprionodon terraenovae 68 and Squalus acanthias 74 . These change in population parameters could provide a compensatory mechanism in small-medium chondrichthyans decreasing the natural and fishing mortality, and likely promoting the recovery of those exploited populations 75 .
In summary, this study shows the feasibility of fishery-independent scientific surveys, which are currently being developed in the Mediterranean, to evaluate the populations trends of by-catch vulnerable species such as demersal chondrichthyans. Our trend models reveal that most of the sharks and batoids currently making up these communities have resisted the impact of fishing on the western Mediterranean over the last two decades. However, this is not the case for two deep water species, which are not resisting fishing exploitation. These results can be explained by the evolution of the trawl fishery (reduction of effort and displacement to deeper waters) over the last few decades, jointly with the great resilience displayed by some of these species, like S. canicula, a small-medium shark which also showed a reduction of its length at first maturity, probably as an evolutionary response to the general overfishing of the demersal populations exploited by bottom trawl fishery since the 1980s in the study area. Our findings, which come from a long period of over-exploitation of these vulnerable species -some of which appeared to be fairly common and have currently almost disappeared or are even considered to be locally extinct -provide a current baseline scenario in order to implement management measures that will strengthen or initiate the recovery of chondrichthyans in the Mediterranean, whose populations play an important role in marine ecosystems, and to develop true adaptive management in the area, making the sustainability of trawl fishery compatible with the conservation of marine ecosystems.  www.nature.com/scientificreports www.nature.com/scientificreports/ Data source. A total of 3158 stations were sampled during the MEDITS bottom trawl surveys, carried out annually during spring and early summer from 1994 along the Mediterranean Iberian Peninsula and the Balearic Islands (Fig. 11). The sampling method is described in MEDITS protocol 76 Supplementary Table S1. The sampling gear was the experimental bottom trawl GOC-73, equipped with a 20 mm mesh size in the cod-end. Hauls were conducted during daylight hours at a towing speed of around 3 knots, with an effective duration of 30' at depths shallower than 200 m and 60' below this depth. The arrival and departure of the net to the bottom, in addition to its horizontal and vertical openings (2.4-3 and 16-22 m, respectively) were measured using a SCANMAR system (http://www.scanmar.no/en/) consisting in a set of sensors attached to the gear, which allow its depth and geometry during the haul to be measured simultaneously.
The catch of each sample was sorted, identified to species level, counted and weighed, and length of individuals was measured. Data on abundance and biomass were standardized to number of individuals per km 2 (referred to as "density") and to kg per km 2 (referred to as "standardized biomass"), respectively. To do so, the bottom surface sampled during the haul was considered, estimated from the horizontal opening of the net, and the distance covered, obtained from a GPS system. Sex and sexual maturity were determined from macroscopic examination of the gonads, following Holden and Raitt's maturity scale 77 , which considers four stages for both females and males: immature, maturing, spawning and post-spawning.
The time series of data on density and standardized biomass spanned from 1994 to 2015 for GSA01 and GSA06, and from 2001 to 2015 for GSA05. Due to changes in the MEDITS sampling protocol over the years, data series on size and sexual maturity were shorter. In GSA01 and GSA06 it comprised the period 2000-2015, whereas available periods in GSA05 were 2001-2015 and 2007-2015 for size and sexual maturity, respectively. Data analysis. Dynamic factor analysis (DFA), a multivariable method applied to relatively short and non stationary time series data to detect common trends 77,78 , was applied to identify common trends between GSAs for diversity, density, standardized biomass and size structure parameters. In DFA, time series data are modelled as a linear combination of underlying common trends and factor loading indicates the strength of the influence of each GSA time series on the corresponding common trend. In order to achieve a proper interpretation of our results, the time series were standardized by subtracting the mean and dividing by the standard deviation 78 . Factor loadings with values above 0.20 (in absolute value) signifies that the corresponding GSA has a high contribution to the common trend detected 77 . Negative factor loading values indicate that the particular GSA has an inverse trend to the one detected as a common trend. Correlation of observation errors was modelled using different covariance matrices: (i) same variance and no covariance (diagonal-equal); (ii) different variances and no covariance (diagonal-unequal); (iii) same variance and covariance (equalvarcov); and (iv) different variances and covariances (unconstrained). The corrected Akaike information criterion (AICc) was used as a measure of goodness of fit, with the best model having the lowest AICc 78,79 . Model implementation was performed using the Multivariate Autoregressive State-Space "MARSS" package 80 in R 3.3.2. www.nature.com/scientificreports www.nature.com/scientificreports/ DFAs were applied to diversity, density and standardized biomass time series in order to search for common trends between GSAs along the study period. In this study, we evaluated the following diversity indices to identify a set of indicators that could provide a good representation of changes in demersal chondrichthyan communities, taking account the different aspects of diversity: mean species richness index (S), the Shannon-Wiener diversity index (H') and evenness (J') that were estimated by year and GSA.
The mean density and standardized biomass of the chondrichthyan groups (sharks, batoids and chimaeras) were calculated by year, GSA and depth stratum, except in stratum A, which was excluded from the analysis due to the low number of sampling stations. These parameters were also estimated by GSA for the most abundant species within their bathymetric distribution range. These species were those representing ≥85% density or standardized biomass in at least two GSAs, or with a percentage occurrence greater than 20% within their depth range of distribution for all GSAs. Since the presence of batoids was practically zero below 500 m depth, stratum E was not considered for this group.
Temporal trends in size structure were also evaluated for the most important species in the three GSAs. To do this, sexes were pooled into different categories according to the total length (TL) of the species, and the mean density by year was calculated. DFAs were applied to identify common trends between GSAs in density values for each length category from 2000 to 2015.
Meanwhile, linear regression analyses were applied to detect temporal trends for the most important species within a particular GSA (i.e. those representing ≥85% density and standardized biomass and with a percentage occurrence within the depth range of distribution ≥20%), in both density and size structure. Dasyatis pastinaca was not included in this analysis because of the possible misidentification with Dasyatis cf. tortonesei in the Balearic Islands 81 .
Linear regression analysis was also used to detect temporal trends of size at first maturity of the most important species. Maturity was assessed for species with at least 50 individuals sampled per year and GSA. To do this, maturation stages were converted into binary data (stages i and ii as immature and stages iii and iv as mature), while length at first maturity (L 50 ; length at which 50% of the individuals are mature) was obtained by fitting a logistic regression to the proportion of mature individuals by length class. This was done for females and males and by year and GSA. The smallest length of mature females per year for each GSA was also estimated. In addition, considering the fact that the condition of fishes has been reported as an important factor affecting growth rate 82 and that an increase in growth rate before the onset of maturation can lead to a decrease age and therefore size and maturity 83 , we also analysed the temporal trend in somatic condition for species with a significant decreasing trend in L 50 . For these species, we calculated the somatic condition for immature individuals (below the length at first maturity reported for these species in the study area 48 ) of both sexes. Individual total weight (TW) and total length (TL) were log-transformed and the linear relationship between logTW and logTL was established by year. Residuals of the differences between the observed and predicted logTW were calculated and standardized, dividing each one by the standard deviation of their predicted values 84 . Lastly, linear regression analysis was applied to detect temporal trends from 2007 to 2015 in each GSA.