|Home | About | Journals | Submit | Contact Us | Français|
Sharks and rays' abundance can decline considerably with fishing. Community changes, however, are more complex because of species interactions, and variable vulnerability and exposure to fishing. We evaluated long-term changes in the elasmobranch community of the Adriatic Sea, a heavily exploited Mediterranean basin where top-predators have been strongly depleted historically, and fishing developed unevenly between the western and eastern side. Combining and standardizing catch data from five trawl surveys from 1948–2005, we estimated abundance trends and explained community changes using life histories, fish-market and effort data, and historical information. We identified a highly depleted elasmobranch community. Since 1948, catch rates have declined by >94% and 11 species ceased to be detected. The exploitation history and spatial gradients in fishing pressure explained most patterns in abundance and diversity, including the absence of strong compensatory increases. Ecological corridors and large-scale protected areas emerged as potential management options for elasmobranch conservation.
Analyses of exploited fish communities in coastal, demersal and pelagic ecosystems have shown that elasmobranch (sharks and rays) diversity and abundance can decline considerably after only short periods of fishing1,2,3. For example, in the Northwest Atlantic, catch rates of 18 large pelagic and coastal shark declined by 49–89% in less than 15 years2. In South Africa, large coastal sharks were reduced by 27% - >99% after 20 years of shark netting programs1, and in Southeast Australia, 20 years of trawling reduced demersal elasmobranch catch rates by >80%3. Industrial bottom trawl fisheries in particular can have strong effects on demersal communities by unselectively catching a wide range of species while destroying complex seafloor habitats4,5.
Aside from declines in population abundance, shifts in community composition can occur because fishes have different intrinsic vulnerabilities to exploitation, are unequally exposed to fishing, and respond to changes in predators and competitors. Any species is characterized by an intrinsic rebound potential (r), a population parameter that combines a set of biological traits (maturity, fecundity, and growth) determining the species productivity and capacity to sustain exploitation6. As a group, elasmobranchs generally have very low r values due to their late sexual maturity, low fecundity and slow growth rate, hence they are unable to sustain even moderate levels of exploitation7,8. Yet even among elasmobranchs, r is variable enough9,7 to sometimes explain observed differences in the species' response to fishing10,11. In other cases, variable exposure to fishing12,13,3 and changes in the abundance of predators and competitors12,14,15,16,17 have been more important predictors of community change.
Predation is an important factor explaining meso-predator population dynamics and community changes in terrestrial18 and aquatic ecosystems19. Among elasmobranchs, large sharks are the main predators of smaller sharks and rays20. While many large sharks have declined over the last decades, meso-predatory elasmobranchs have increased in exploited communities of the Atlantic21,17,15,12, Pacific22,3 and Indian Ocean23,24. However, the magnitude and trajectory of population increases are not always predictable because prey benefit from both reduced natural mortality and a reduced risk of predation25,18. The latter is a more pervasive but not easily quantifiable behaviorally-mediated effect that influences prey distribution, habitat and food choices, and overall fitness25. Moreover, while sharing common predators, sympatric species compete for limiting resources, which influence their abundance. When fishing or other human impacts selectively deplete or removes species26, releases from competition are often predicted and were suggested to explain unexpected increases in elasmobranchs in several areas, of the world21,16,17,27. Yet in some of these cases, alternative explanations, such as shifts in spatial distribution, also played a role28.
Explaining community changes in large marine ecosystems by identifying the occurrence and relative contribution of direct and indirect effects of fishing, or other human impacts, is challenging because it requires performing controlled experiments impractical at large scales. In these cases, adopting an observational approach across gradients of natural or human-induced perturbations is an efficient alternative19. We used such an approach to study the long-term changes in the elasmobranch community of the Adriatic Sea, whose long history of human-induced changes, as well as the spatial and temporal contrasts of perturbations, offered an ideal case study. The Adriatic Sea is a semi-enclosed Mediterranean basin that has been exploited for thousands of years, and where large marine predators (sharks, pinnipeds, and cetaceans) have declined dramatically over the past two centuries29. The Adriatic's broad continental shelf and accessible fishing grounds allowed the development of large fisheries for shellfish and groundfish30. Yet fisheries developed unevenly between the western and eastern sides of the basin. While Italian waters were exposed to extremely high exploitation pressure from high-capacity fishing fleets, former Yugoslavian sectors sustained a much lighter fishing exploitation until recently (see Supplementary Methods for a brief fishing and ecosystem depletion history).
To examine spatial and temporal changes in the elasmobranch community, we used data from five different scientific trawl surveys carried out in the Adriatic since 1948 (Fig. 1, Table 1). First, we extracted species-specific trends of catch rates to estimate short-term community changes within surveys. Then we estimated long-term community changes comparing catches across surveys. Finally, we used life-history characteristics, environmental factors (e.g. sediment composition, temperature), fishing effort data, and historical fishing information to explain the observed trajectories of population and community changes. It was our aim to identify the main drivers of change, whether the decline of large predatory sharks triggered meso-predator releases, and if there was a change of the elasmobranch community composition suggesting competitive interactions.
In total, 2575 tows carried out over six decades across the Adriatic Sea detected 33 small, demersal, meso-predatory elasmobranch species (average trophic level: 3.9, SD: 0.1231), including 12 sharks, 20 rays and one chimaera (included for their evolutionary and ecological similarity, Table 2). Of these, 11 species ceased to be detected during the period of observation (no more occurrences after the year 2000, Table 2), while 6, mostly deep-water species and small skates, were only recently detected by the MEDITS surveys which expanded to greater sampling depths (Table 2).
Across all five surveys, species-specific frequency distributions were very skewed (Fig. 2), with few dominant species and many occurring only sporadically. The Hvar survey detected the highest diversity (23 species, Shannon Index [SI]: 3.39, Fig. 2a), dominated by small-spotted cat sharks and thornback skates with unstandardized densities of 426.8 and 76.8 individuals/km2, respectively, and frequencies of occurrence (FO) of 0.76 and 0.71. The other 21 species were caught in < 21% of the tows with densities below 11 ind./km2. Over time, moving to the most recent survey, richness and abundance decreased toward more flattened and truncated distributions. In MEDITS (SI: 1.96), the small-spotted catshark was still the most abundant species, but with a density of 62.1 ind./km2 and FO of 0.20 (Fig. 2), a 6.8- and 3.8-fold decline, respectively, compared to Hvar. Overall, 21 out of 27 species had FO < 0.021 and densities < 4 ind./km2. The high elasmobranch abundance and diversity characterizing the central Adriatic during the Hvar survey in 1948–49 disappeared (Fig. 3a and b). Yet, species richness and abundance were higher in the eastern coastal areas than elsewhere (Fig. 3a and b). Elasmobranch abundance in Croatia was almost one order of magnitude higher than in Italy (Supplementary Fig. S1), where sharks and rays were largely absent except for a relatively high-density zone in the upper Adriatic (above the 50 m isobath, Fig. 3a MEDITS) mainly composed of spurdogs (Supplementary Fig. S2), smooth-hounds, and eagle rays.
Standardized catches generally confirmed the above patterns. Although the fitted models had considerable selection uncertainty - in each survey, more than half of the species had a selected best model with <10% chance (Akaike weight <0.1) of being the most plausible among a 95% confidence set of models (see methods, Supplementary Tables S2–S6) – spatial covariates (depth, latitude and longitude) resulted as the best predictors of species abundance (Tables 3 and S2–S6). In the largest surveys available (Hvar and MEDITS), these predictors pointed to an increase in catch rates from the Italian to the Croatian coasts, and for most species from deep to shallow waters (i.e. from offshore to inshore, Supplementary Fig. S3). Temporal covariates (year, time of the year) were less important than spatial variables, but many species had significant short-term temporal trends within surveys (Fig. 4).
Analyses of catches in the Jukic survey (1963–1971), located in the central eastern Adriatic, identified nine species with reliable short-term trends in abundance, and all except the thornback skate showed an increase in standardized catches (three species were significant, Fig. 4a). In comparison, between 1994 and 2005 across the Adriatic (MEDITS surveys), 16 species had reliable estimates of population change. Of these, nine species, mainly sharks, showed declines (3 statistically significant, Fig. 4b), while increases were mostly shown by meso-pelagic rays and small skates; yet these trends were not significant except that for the eagle ray, which increased by 3.74 times (Confidence Interval, CI: 1.06, 13.45) in 11 years.
Over longer periods of time, changes in abundance were larger and more significant. Comparing Hvar and MEDITS surveys showed that elasmobranchs declined by 94.5% over 57 years (Fig. 5a and b). Sharks declined more than rays (−95.6% vs. −87.7%) with small-spotted catsharks (−96.2%, CI: −97.8%, −93.5%) driving most of the patterns. Rays shifted in species composition. In particular, the thornback skate, the most abundant ray in the 1940s, recorded the steepest decline (−97.2%, CI: −98.4%, −95%), whereas brown skates increased by 2.36 times (CI: 1.05, 5.3) becoming the most abundant skates (Fig. 2e). Significant long-term increases were also detected for eagle rays (111 times, CI: 17.05, 735), marbled torpedoes (75.9 times, CI: 5.07, 1135.64) and spurdogs (3.1 times, CI: 1.05, 9.27) (Fig. 5a and b).
In the coastal area off the eastern Adriatic (Zupanovic area, Fig. 1), between 1957 and 2005 elasmobranchs increased significantly by 2.12 times (CI: 1.59, 2.83). Rays increased faster than sharks (3.79 vs. 1.88 times) with brown skate and common eagle rays recording fold increases of 17.7 (CI: 9.9, 31.7) and 19.73 (CI: 7.08, 55.03) respectively. As for sharks, the patterns were largely driven by the almost 2-fold increase of small-spotted catsharks (1.92, CI: 1.1, 3.4). The smooth-hound recorded the largest increase (21.13, CI: 8.23, 54.25), and the thornback skate was the only elasmobranch showing a significant decline (−39.8%, CI: −63.14%, −1.77%) (Fig. 5c).
In the Jukic area, comparing catch rates in the period 1948–1971 (Hvar-Jukic comparison), elasmobranchs declined by 61% in 23 years. These were mainly sharks (−61%), while rays recorded a moderate and non-significant decline (−38%), although brown skates suggested some increase, and the common eagle ray increased by 128 times (CI: 46, 358). The spurdog was the only shark increasing significantly (4.43-fold, CI: 1.43–13.72), while 4 of the other 5 species we could model recorded declines between −63% (small spotted cat-shark, CI: −83.6%, −16.6%) and −96.8% (the smooth-hound, CI:−98.7, −92%). In contrast, in the period 1963–2005 (Jukic-MEDITS comparison), sharks and rays as a group declined comparably and significantly by more than 90%. Spurdogs reversed their earlier increase with a decline of −87.6% (CI: −95.2%, −68%) and the increase of common eagle rays became non significant. This time, smooth-hounds recorded the only significant increase (356 times, CI: 6.8, 18534; Fig. 5e).
All these results were robust to mis-specification and uncertainty of trawl performance (in terms of swept area), to sampling and estimation error of sediment composition, and to the surveys' sample sizes (see Supplementary Methods).
Our projected distribution of current fishing intensity (Fig. 2f) reflected the patterns in abundance and distribution detected in MEDITS (except for the upper Adriatic Sea, Fig. 3), and was consistent with the spatial parameter estimates of the standardized catches (Supplementary Fig. S3). The eastern Adriatic (mainly fished by Croatia) has a much lower level of fishing pressure than the western side (exploited by Italy). Italy records about twice the amount of otter trawlers (1541) than Croatia (855) with an average trawler having about 2.25 times the Horse Power (199 HP, sd = 149) of an average Croatian one (88.5 HP)32. We predicted a higher fishing intensity all along the northwestern Adriatic, especially between Fano and Pescara, and around Chioggia (Fig. 2f). While in Croatia, heavily fished areas would be the southeastern part of Istria, and above the Dalmatian channels between Sibenik and Split. Note that the distribution of fishing effort was coarser in Croatia because our source data was aggregated by major fishing districts, and we assumed that the allocation of fishing boats per port was proportional to the port population density (Supplementary Materials). Dividing the number of boats fishing in these different Croatian sectors by the sectors' trawable surface (i.e. excluding 3 nautical miles from the shores) showed that most of the Croatian channel region would be free from trawling, which instead would take place mostly in the coastal areas (between 3 and 6 nautical miles) of Istria and the outer channel areas (Supplementary Fig. S4). If the number of otter trawlers longer than 18 meters is a good index of offshore fishing intensity (between 3 and 40 nautical miles), fishing exploitation is an order of magnitude higher in Italian than Croatian waters (0.005 vs. 0.0003 boats per square kilometer, Supplementary Fig. S4).
The intrinsic vulnerability of species did not explain the variability of the observed rates of change. We found no significant relationships between r and βys estimated from MEDITS (slope 0.996, p-value: 0.504, R2: 0.035) Jukic (slope = 3.499, p-value 0.152, R2: 0.27), historical comparisons in the Hvar area (slope = 14.020, p-value: 0.51, R2: 0.039), Jukic area between Hvar and Jukic (slope = 5.91, p-value: 0.75, R2: 0.013), Jukic area between Jukic and MEDITS (slope = 14.53, p-value = 0.45, R2: 0.072) and in the Zupanovic area (slope = 23.61, p-value: 0.18, R2: 0.21).
Understanding long-term changes in exploited fish communities requires the consideration of several factors, including the intrinsic vulnerability of various species to exploitation, changes in biological interactions (e.g. predation and competition), different exposure to fishing (e.g. catchability, availability, and commercial value) and susceptibility to other stressors such as habitat degradation and pollution. All of these factors can alter the species-specific response to exploitation, generating complex community changes over time and space1. Here, we attempted to understand long-term changes in the elasmobranch community in the Adriatic Sea and the underlying drivers.
By analyzing trawl surveys carried out in the area over the last six decades, we detected a structurally depleted elasmobranch community. A total of 2575 tows were not sufficient to detect at least 25 elasmobranch species that were recorded in the area before 1948 (Supplementary Table S7). These included 13 great sharks (>2 m in length), which may not have been caught due to their lower trawl catchability compared to demersal species, but also 8 bottom associated large meso-predatory elasmobranchs (e.g. angelsharks, guitarfishes, large rays, bramble and gulper sharks, Supplementary Table S7), some of which were the target of dedicated fisheries33,34. Most of the 33 species detected by the trawl surveys strongly declined over time, and 11 disappeared, updating the list of possible extinctions in the Adriatic Sea14,1 to 22 species. Overall, sharks declined stronger than rays (−95.6% vs. −87.7%), and more shark than ray species recorded significant declines (4 vs. 2 species respectively, Hvar-MEDITS comparison Fig. 5b). Similar strong depletions of elasmobranch abundance and diversity were recorded in exploited fishing grounds of the North Sea27, southeastern Australia3, the Tyrrhenian Sea35, and the Gulf of Lions36. However, none of these studies covered trends across six decades with standardized fishery independent data.
The Mediterranean Sea in general and the Adriatic Sea in particular, have experienced a long history of human impacts, including exploitation, pollution, and habitat degradation, resulting in severe population declines of marine species37,29. With our data starting in 1948, we only captured a final stage of this depletion history, with many large sharks already depleted or absent. However, we found a higher abundance and diversity of elasmobranchs in the eastern Adriatic (Fig. 3a and b) which reflected the less intense historical and recent fishing pressure in Croatian compared to Italian waters (Fig. 2f), especially in coastal Croatian waters were trawling has been restricted, discouraged or hardly practicable38 (also see Supplementary Fig. S4). This good correspondence between the spatial gradient of fishing pressure and the observed patterns of elasmobranch abundance and diversity suggests that fishing has been a strong driver of the observed changes. Other human impacts, such as pollution and habitat degradation, were likely also stronger on the more highly populated Italian coast and may have contributed to the observed declines.
Differences in the intrinsic vulnerabilities among species (reflected by their rebound potential r) were not sufficient to explain the observed species-specific rates of change. They were suggestive of a positive but not significant relationship between vulnerability and rate of decline. The surveys detected a limited and biased sample of species whose life histories did not represent the full spectrum characterizing the native Adriatic elasmobranch fauna. In 1948 (when our trawl surveys began), coastal predator communities were already depleted37,29, and demersal elasmobranchs were nearly absent on commonly trawled grounds, even on the eastern side39. Although offshore grounds remained almost unexploited until the end of WWII (Supplementary Materials), the surveyed fish community was characterized by small, productive elasmobranchs. Larger, less resilient meso-predatory species (common skates, tope sharks, angel sharks, and others), which used to be common or seasonally abundant throughout the basin in the 19th and early 20th century33,40,41,42, were already scarce or under detectable levels (present in <5% of tows, Table 2, Fig. 2a) likely due to decades of directed and incidental coastal fishing33,34. Some species of low commercial value (e.g. torpedoes and eagle rays)43 might have benefited from a lower exploitation rate, but few elasmobranchs are discarded in the Adriatic if reaching a marketable size40,44,43.
The observed species-specific trends were likely not independent of changing interspecific interactions such as predator or competitor releases. Although compensatory changes in population abundances in the Adriatic were less visible than in other marine regions such as in the north west Atlantic and Gulf of Mexico12,15, some increases occurred in the least exploited eastern sectors in historical surveys (e.g. Jukic area Fig. 4a), and over long periods (e.g Zupanovic area Fig. 5). In general, competition and predation releases are more evident when not confounded by high levels of fishing mortality12. Thus, when we compared historical and more recent surveys off central Croatia (Jukic area, Fig. 5d and e), an area where offshore fishing pressure has only recently increased38, earlier population increases were later reversed or buffered such as for spurdogs and common eagle rays. Similarly, long-term patterns of change detected over the Hvar area, were different in the less exploited Croatian channels during a comparable time span. There, only thornback skate recorded a significant albeit moderate decline (−39.3% over 48 years, CI: −63.1%, −1.8%) while most of the other species increased (Zupanovic area, Fig. 5 b and c). Thus, compensatory changes might persist in relatively unexploited areas while in other parts were eventually overruled by the effects of exploitation.
Stronger predator or competitor releases might have occurred historically in the 19th and early 20th century, when large predatory sharks (e.g. hammerhead, mako, porbeagle and white sharks) were in decline37,45, but fishing had not yet expanded to industrial levels (Supplementary Methods). At that time, angel sharks, spurdogs, smooth-hounds, and skates were so abundant to sustain targeted fisheries33,34, and the main Adriatic fish markets recorded increases in elasmobranch landings44,46. Such increases were particularly evident after periods of reduced fishing operations (e.g. during the two World Wars) and characterized by exceptional but not biologically plausible catches. For example, in the years immediately following the wars, D'Ancona (1949) described foot-seine hauls of 70–160 smooth-hounds averaging 10 kg a piece in Chioggia, despite the fact that it takes about 16 years for that species to reach such a size (Supplementary Methods). This may suggest compensatory population increases of residential sharks, but also possible immigration of individuals from neighbouring, less-exploited areas. Similar patterns were observed in the North Sea47, and on Georges Bank (NW Atlantic) where the rapid increase of winter skate, previously attributed to competition releases, were related to environmentally driven immigration of individuals from the northern Scotian Shelf28.
Differences in species mobility, and thus exposure to exploitation are congruent with the observed patterns of resilience and depletion in the Adriatic. The high densities of bentho-pelagic elasmobranchs (spurdogs, smooth-hounds, and eagle rays) in the heavily exploited upper Adriatic might partly highlight the simplistic nature of our fishing-effort model, but also a possible spillover of mobile species from the underexploited Istria to the heavily fished but extremely productive northern Italian sectors. Specific oceanographic conditions, the influence of the Po River, and strong seasonal fluctuations of oceanographic regimes30 make the upper Adriatic an extremely productive area characterized by a great abundance of benthic (molluscs and crustaceans) and pelagic primary consumers (sardine and anchovies)30, which attract elasmobranchs. Devoid of pelagic larval stages, elasmobranchs disperse relying on intrinsic mobility and available ecological corridors48. The species that were dominant in the upper Adriatic can undertake transoceanic and long seasonal migrations49,48,50, and have home ranges much larger than the northern Adriatic. Hence they may be little affected by localized stressors because they are able to replenish from areas even outside the Adriatic, possibly using relatively sheltered gateways in Croatian coastal waters. Conversely, elasmobranchs with limited mobility such as catsharks and other skates (home ranges <50 km)47,49 persisted in the less-exploited eastern Adriatic and almost completely disappeared from the heavily exploited Italian waters.
In summary, the long history of human impacts severely depleted the Adriatic predator and meso-predator community, leaving the ecosystem in a structurally altered state. The elasmobranch community went through a sequence of population depletions (from great sharks in the 19th and 20th century37, to large meso-predatory elasmobranchs in the early 20th century to smaller more productive species in the later 20th) and shifts in species composition. More recently, this culminated in the near disappearance of elasmobranchs in the most exploited areas of the Adriatic Sea. Cases of ecological compensation due to release from predation or competition were not as evident as in other marine regions and may have pre-dated our observation window. Yet the strong gradient of fishing intensity from the Italian to the Croatian side allowed the persistence of a more abundant and diverse elasmobranch community in the eastern Adriatic, which may fuel a spillover of mobile species toward more exploited western areas. Similar contrasts of decline and resilience of exploited fish populations were observed in large systems of marine protected areas established on Georges Bank and the Great Barrier Reef51,52,53.
In January 2008, Croatia established a 23,870 km2 Ecological and Fisheries Protection Zone (EFPZ) in its national and international waters54. However, after 3 months, the EFPZ was dismissed because of harsh opposition of bordering countries. Our results suggest that planned and managed internationally, and concerned with a reduction of fishing pressure in western areas and control of developing fisheries in the east, a similar initiative could protect abundance and diversity of Croatian resources from further depletion. Furthermore, it could promote elasmobranch recovery in the overall Adriatic by profiting from the remaining strongholds of elasmobranch abundance and diversity in the east, otherwise scarcely protected by current marine protected areas55.
Our dataset comprised 2575 trawl tows, from three surveys identified as Hvar, GRUND, and MEDITS, covering large portions of the basin; and two surveys, we called Jukic and Zupanovic, more locally confined to Croatian waters (Fig. 1, Table 1, Supplementary Methods and Table S1).
We initially analyzed each survey individually to avoid dealing with differences in sampling framework and gears, potentially influencing the estimated indices of abundance. We identified the available environmental and methodological covariates explaining most of the catch variability, and used these variables to standardize the catches. We estimated short-term trends in catch rates from surveys having more than three years of observations, and estimated long-term changes by combining historical and recent surveys in overlapping spatial domains. Finally, we attempted to explain the observed spatial and temporal patterns of the catches with data on fishing effort and species life-histories.
Standardization is paramount when combining independent experimental observations (such as this set of surveys) differing in sampling framework, gear specification and survey's spatial domain (area and depth range). This process avoids inferring false population trends from observed catch trajectories. Also taking into account the inherent density, detection probability, and spatial distribution of a species is necessary to prevent estimating biased indices of abundance. The low population densities of elasmobranchs and their patchy distribution generate catch distributions highly skewed, successfully described with a negative binomial distribution2,12,37.
GLMs flexibly accommodate unbalanced sampling schemes and non-normal data56. Hence using this framework, we assumed that the chance of obtaining n individuals of a species in any towi followed a negative binomial distribution with mean µi and variance ; k is a dispersion parameter estimated from the data. We modeled log(µi) as a linear function of a number of covariates characterizing each tow, log(µi) = α + XB + log(Ai), where α is the intercept, X the matrix of covariates, B the vector of their relative parameters, and Ai is the swept area, treated as an offset. Ai was usually measured by multiplying the trawl net horizontal opening (ho) by the distance towed (details in Supplementary Methods).
For each survey, from a set of available covariates (Supplementary Methods) we selected a subset (Table 3) to construct an initial model from which we proceeded with model selection. For avoiding computational problems, arising from fitting variables with different numerical scales57, continuous variables were standardized by using unit normal scaling, , where is the mean of variable j, and sj is its standard deviation. To avoid collinearity, correlated variables were discarded or combined into predictors expressing unrelated information.
We needed to identify important predictors, answer ecological questions, and also avoid the identification of spurious relationships resulting from over-fitting a full-saturated model (having 2n parameters, for n variables, including all interactions of main effects) to small datasets. Therefore we included only main effects and quadratic forms of some continuous variables (Table 3). Quadratic functions of depth, temperature, latitude and longitude were included to capture the tendency of animals to aggregate around optimal values of environmental variables. A quadratic function of year in MEDITS, the survey with the longest time span, was also included to model possible compensatory responses of species to changing exploitation regimes, competition or predation. Higher order temporal polynomials were not attempted, having limited temporal observation windows on animals characterized by slow population dynamics. To capture the habit of many elasmobranchs to undertake in seasonal migrations or shift in distribution58, we included a periodic function of survey date, , where Jdi is the ordinal day for tow i (referred to the earliest date of all the surveys), which eventually allowed us to compare different surveys for tows performed at different times of the year. Substrate composition was included to account for the strong reliance of demersal species on seabed features59.
From the initial model, for each species, we selected a reasonable preliminary model structure by backward elimination of covariates in SAS 9.1 (Supplementary Methods). This preliminary model was used to estimate the negative binomial dispersion parameters k to use in a second model selection stage (performed in R 2.15.060), when most variable combinations of the initial model (number is given in Table 3) were refitted by providing the estimated k's as initial values. For each model mi (variable combination) we calculated the Akaike Information Criterion (AIC), AIC differences (Δi = AICi – AICmin), Akaike weights , where R is the number of models fitted; and evidence ratios (wmax/wi)61. We then selected the best model corresponding to the minimum AIC value, and the set of models with AIC differences < log(1/8) * 2 to identify a 95% confidence set (CS) of other plausible models61. We used this CS to calculate the most important variables affecting the variability of the species, and the overall importance of these variables in explaining the catches of most species. To calculate the importance (w+) of a variable (xj), we summed the Akaike weights of all the models containing the variable among the CS, , where Ij(mi) = 1 if variable xj is in model mi, 0 otherwise; R* is the number of models in the CS. To evaluate the overall importance of the variables across species, we first ranked the variable importance for each species from 1 (the most important, highest w+) to n (n ≤ the maximum number of covariates used in a given survey). Then we averaged the ranks across all species.
In Jukic and MEDITS, we estimated a species-specific instantaneous rate of change over time βy by including year among the set of standardizing covariates regardless of whether it was included in the best model. We profiled the best-model likelihood functions for a range of fixed βys (−0.9, 0.9, i.e from >99% declines to about 20,000 fold increases) and selected the value corresponding to the maximum likelihood. Obtaining within-survey trends was important to test the species' response to current exploitation levels and uncover possible compensatory community changes due to competition.
We combined spatially overlapping surveys to estimate long-term changes in CPUEs for all elasmobranchs combined as well as for all sharks (including rabbit fishes), rays, and individual species. We compared the tows performed in overlapping survey domains (i.e. using survey areas as Venn diagrams and comparing tows within intersections) to avoid unbalanced spatial comparisons. Of all the possible intersections, four were instrumental to evaluate changes in the species abundance and composition under different regimes of fishing: 1) in the Hvar domain (Fig. 1), we compared Hvar and MEDITS tows excluding those in MEDITS deeper than the deepest tow in Hvar (1948–2005); 2) in the Jukic domain, we compared Hvar and Jukic tows (1948–1971), and 3) Jukic and MEDITS tows (1963–2005); 4) in the Zupanovic domain, we compared Zupanovic and MEDITS tows (1957–2005). For all comparisons, we used a common glm structure with negative binomial error distribution and log link, where the expected catch per tow (µi) was a function of a categorical variable (SV) indicating whether towi was performed in the old or recent survey, the most important variables identified in the previous analytical stage, and time of the year to account for seasonal differences across surveys, η = log(µi) = α + βsvSVi + βs1S1i + βs2S2i + βnNi + βeEi + βdDi + log(Ai); covariates are defined in Table 3. The initial dispersion parameters k's were those estimated in the standardization stage for the MEDITS surveys.
If observed patterns of change in the exploited community are independent of interspecific interaction, habitat degradation or other factors affecting the population dynamics of the exploited species, we should expect a correspondence between the species' intrinsic vulnerabilities and their rates of change. Therefore we fit weighted regressions (for each short- and long-term trend analyses) between the species-specific changes in catch rates (βy and βsv), and the species' r estimated following Smith et al.9. As weights, we used the inverse of the parameter estimates' variance. Life-history parameters necessary to estimate r were collected from multiple bibliographic and electronic sources listed in the Supplementary Methods.
Similarly, as fishing intensity is spatially heterogeneous across the area, we should predict to have higher elasmobranch density in less exploited areas and vice versa. Therefore, we visually compared the spatial distribution of CPUEs (catch per unit effort) with the distribution of a predicted index of current trawl fishing intensity projected over the Adriatic. By using data on number, horse power (HP), and gross tonnage (GT) of boats fishing along the Adriatic coasts, we assumed that the intensity of fishing per unit of sea bed surface (at any point in the Adriatic Sea) was proportional to the cumulative effort exerted by all Adriatic fishing fleets, whose individual contribution depended on the distance between any point and the fleet harboring port, fleet size, and point bathymetry. Using empirical data of effort distribution of a sample of four Italian trawl fleets we developed a working model predicting fishing intensity (daily units of HP per unit area) from bathymetry, distance from port and fleet size (data sources and analysis details are found in the Supplementary Methods). This analysis only provided an approximate spatial picture of fishing effort. We omitted local scale predictors of fishing intensity such as avoidance of fisheries' competition, resource abundance, compliance with spatial and seasonal management, and geopolitical restrictions. Nevertheless, for our purposes of distinguishing general patterns of fishing effort between the eastern and western, northern and southern parts of the Adriatic, our analysis should provide a sufficient approximation.
F.F., G.C.O. and H.K.L. conceived and designed the study; F.F., G.C.O. and C.J.J. collected and analyzed the data; F.F., G.C.O., A.A.R. and H.K.L. wrote the manuscript.
Supplementary Methods and Results
We would like to thank P. Casale for providing data on empirical effort distribution, the Italian Ministry of Agriculture and Forestry for granting permission to use MEDITS data, C. Piccinetti for providing MEDITS data and W. Blanchard and E. Susko for statistical advice, F. Serena for sharing data on trawl performances and discussions, F. Micheli, C. Mora, N. Dulvy, G. De Leo, Chris Fehily and Lauren Daly for comments and discussions. Financial support was provided by the Lenfest Ocean Program and the Oak Foundation.