Analysis of Long-Term Changes in a Mediterranean Marine Ecosystem Based on Fishery Landings

Introduction

Long-term analyses of the interactions between human society and the oceans are necessary for understanding the processes that brought the marine ecosystems as we see today, avoiding the so-called “shifting the baseline syndrome” and understanding the magnitude and causes of change (Pauly, 1995; Jackson et al., 2001). In this framework, marine historical ecology (MHE) can bring a significant contribution to present-day management of marine ecosystems, both for conservation and for sustainable exploitation of resources (Engelhard et al., 2015). However, often the analyses of historical changes in marine communities cannot be based on results of structured and standardized monitoring programs, since in most of the cases those were set in the very last decades. The need to bridge the gap between request of knowledge on past status of ecosystems and the available monitoring data, triggered the uses of different approaches to extract information from data coming from other sources, including paleontological, archeological, and historical sources, as well as the use of landings (e.g., Rosenberg et al., 2005; Sàenz-Arroyo et al., 2005; Lotze et al., 2006; Fortibuoni et al., 2010, 2016; Van Beveren et al., 2016).

In this context, detailed and disaggregated fishery statistics represent an important source of information that can be used as proxies to evaluate long-term changes in marine fisheries and communities. Changes in the composition of landings, evaluated through opportune weighting metrics (indicators), showed to reflect changes in the structure of underlying fish communities due to anthropogenic impacts and environmental changes (e.g., Caddy, 1993; Pauly et al., 1998; de Leiva Moreno et al., 2000; Pinnegar et al., 2002; Libralato et al., 2004; Pauly and Watson, 2005; Baeta et al., 2009; Munyandorero and Guenter, 2010; Kleisner et al., 2013). Moreover, catch statistics are recognized to be linked to fishing and environmental pressures and respond selectively to management action (Coll et al., 2016).

Nevertheless, the intrinsic limitations of fishery-dependent data—such as landings—include the lack of standardization, the dependence from fishing fleet activity features as well as market preferences of products that usually change across time and space. All this imposes caution in deriving marine population densities directly from catch statistics (e.g., Essington et al., 2006; Pauly et al., 2013).

The capability to connect modification in landings composition to changes in the community at sea is more robust when local disaggregated landings result from multi-target and multi-gear fisheries (i.e., several distinct métiers), and when changes in fishing activities (e.g., introduction of new technologies, shift from one fishing gear to another, shift in fishing grounds, fishing capacity as number and tonnage of boats) are traceable. In this context an important aspect to be considered is the trend in indices over time, rather than the absolute values they assume, being reference points or limit values for many indicators not yet been established (Shin et al., 2010).

The Northern Adriatic Sea (Mediterranean Sea) represents a valuable case study for MHE, due to the long history of exploitation (Botter et al., 2006; Fortibuoni, 2010), human-induced changes (Lotze et al., 2011), as well as documented long-term modifications of the physical and chemical oceanographic characteristics of the basin due to anthropogenic impact (Mozetič et al., 2009; Solidoro et al., 2009) and temperature change (Russo et al., 2002).

In this study, a long-term time-series (1945–2014) of landings disaggregated by species (Mazzoldi et al., 2014) for the Northern Adriatic Sea was analyzed to detect changes in total yields and variations in landings composition by functional groups over time. A suite of ecological indicators was applied to landings data to integrate responses to multiple stressors (Fu et al., 2015; Coll et al., 2016). They include trophodynamic indicators (e.g., mean trophic level, primary production required), climatic indicators (mean temperature of the catch), and other indicators, such as elasmobranchs-bony fish ratio. These indicators were compared with independent data describing main ecosystem drivers, in order to corroborate findings.

Our analysis falls within the general need of taking into account ecological processes when considering long-term changes in fishery resources where fishery-independent data are lacking, by testing several ecological indicators, and their responsiveness to fishery and environmental drivers. The approach is suitable to be applied for the purposes of the Ecosystem Approach to Fishery Management (EAFM) and contributes to establishing historical baselines to be used to compare current and future ecosystem status, for instance in the context of the Marine Strategy Framework Directive (MSFD) implementation.

Materials and Methods

Area of Study

The Northern Adriatic Sea is the shallowest (average depth of 33.5 m) and northern-most area of the Adriatic and Mediterranean Seas (Figure 1). This area is characterized by strong riverine outflows from the Po River (the largest Italian river) that is a primary source of nutrients and organic matter to the basin (Giani et al., 2012). Water circulation is dominated by two counter-clockwise gyres, which confine a large part of the nutrient-enriched riverine inputs along western coastal regions (Zavatarelli et al., 1998). In fact, the western coastal waters were long considered as ones of the most productive of the Mediterranean Sea (Hopkins et al., 1999), and occasionally local hypoxia/anoxia events have been reported especially close to the Po River mouth (Giani et al., 2012). The reduction in phosphorus loads in Italian rivers in the 1980s triggered reversal in the eutrophication trend and was indicated as the start of a (cultural) oligotrophication process for the basin (Mozetič et al., 2009; Solidoro et al., 2009).

Figure 1. The area of study (Northern Adriatic Sea). Commercial fisheries landings derive from the area shown in the figure, Slovenian and Croatian waters excluded.

The seabed is characterized by muddy and sandy bottoms, with the presence of few rocky outcrops. Due to the presence of a wide flat trawlable platform coupled with the high productivity, the Northern Adriatic Sea is the Italian basin with the highest fishery pressure, and one of the most exploited areas in the Mediterranean Sea (AdriaMed, 2004).

History of Fisheries in the Northern Adriatic Sea

The biological resources of the study area have been intensively exploited since centuries, and the port of Chioggia (Figure 1) hosts the most important fishing fleet of the Adriatic Sea and one of the most important of the entire Mediterranean basin (Botter et al., 2006; Fortibuoni, 2010; Mion et al., 2015).

Before and immediately after World War II (WWII) the Chioggia fishing fleet adopted a wide variety of artisanal gears, both passive and active and the industrialization of fisheries in the Adriatic Sea gradually became consolidated starting from the 1950s, having begun in the period between the two world wars with the first experiments with engines (Fortibuoni, 2010). The period following WWII was characterized by marked changes in fishing equipment and technologies. There was a great increase in both demand for fish products and technical innovations, which enabled fleets to expand considerably. Consequently, there was a continuous and substantial increase in the fishing capacity (Cataudella and Spagnolo, 2011).

The use of progressively larger ships and engines allowed fishing areas to be expanded as well as larger and heavier gears to be used. As an example, the bottom otter-trawl was previously towed by pair sailing boats (Botter et al., 2006), and after WWII it was increasingly adopted by single mechanized vessels. Most traditional fishing gears were abandoned (e.g., beach seines and longlines, commonly used throughout the nineteenth century) and new and more efficient ones were introduced. However, the geographical characteristics of the Northern Adriatic Sea (shallow and semi-closed basin) constrained the expansion of fisheries, thus fishery grounds exploited by Chioggia's fleets can be considered almost stable throughout the period analyzed (Mazzoldi et al., 2014).

The first of the major innovations in fishing equipment introduced was the “saccaleva” surrounding net aided with a light source for attracting fish. It became widespread in the 1940s in Chioggia and gradually replaced all other methods for fishing small pelagic fish, such as the traditional “menaide” drift net (Fortibuoni, 2010). Successively, the “saccaleva” was substituted in the late 1960s by the more efficient “volante” (mid-water pelagic trawl, towed by paired vessels) that since the 1970s is the main gear used to catch pelagic species (Cingolani et al., 1996). In the mid-1950s the “rapido” trawl (a sort of beam trawl rigged with 10 cm long iron teeth; Pranovi et al., 2001) was introduced targeting flatfish and shellfish (such as the Mediterranean scallop Pecten jacobaeus), but also demersal resources, as spottail mantis shrimp (Squilla mantis) and common cuttlefish (Sepia officinalis). In the early 1970s, the first hydraulic dredge came into use, substituting the traditional hand-maneuvered gears to harvest different species of clams (Romanelli et al., 2009). In relation to this, the clam fishery quickly became one of the most valuable.

Furthermore, in the mid-1980s, the introduction of LORAN (Long Range Navigation) and subsequently the video plotter and GPS (Global Positioning System) greatly improved navigation precision, allowing the exploitation of areas that were previously inaccessible because of their proximity to unsuitable trawling sites, as the presence of rocky outcrops.

As regard fisheries management, a comprehensive scheme in Italy was initiated with the Law 41/1982, establishing that all professional fishing vessels had to possess a license reporting the characteristics of the vessel (e.g., GT), limitations of fishing areas, gear use and spatial licensing (Piroddi et al., 2015). The Italian fisheries management system is actually based on fishing effort/capacity regulation systems, and technical measures. No quotas or TACs (total allowable catch) have been established, except for Atlantic bluefin tuna (Thunnus thynnus) and Striped venus clam (Chamelea gallina). From 1983, the national fishing fleet was subject to reduction constraints in relation to two reference parameters, fleet tonnage and engine power. A further incentive toward fleet capacity/effort reduction was provided by European Structural Funds that financed the voluntary removal of vessels. Moreover, in recent years there has been a voluntary departure from the sector due to the general old age of the fishing fleet and the fisheries crisis (Cataudella and Spagnolo, 2011). However, as regards the Chioggia's fishing fleet, fishing capacity (i.e., total GT) has increased until recent years (Barausse et al., 2011).

Data

Landings Dataset

Official landings data (1945–2014) from the Chioggia's wholesale fish-market were retrieved from the Clodia database (Clodia database, 2015). Landings of a wide variety of benthic and pelagic species represent the aggregated commercialized quantities caught by the highly diversified fishing gears employed by fishermen of Chioggia. Data do not include any estimate of the discard, and landings disaggregated by gear are not available. Landings refer only to fish and seafood caught by local fishermen from the Chioggia's fleet that operates in the Adriatic Sea. The database was validated by Mazzoldi et al. (2014) to show that landings composition provides reliable indication of fish abundance.

Information on the habitat (pelagic/demersal) and maximum length (L_max) for each species were taken from FishBase (Froese and Pauly, 2016). Moreover, the thermal preference (median temperature preference, T) was assigned according to Cheung et al. (2013). The trophic level of species (TL), specific for the Northern Adriatic Sea, was obtained from Fortibuoni et al. (2013), from other literature and estimates (e.g., Stergiou and Karpouzi, 2002) or retrieved from FishBase in the case of fish, and SeaLifeBase (Palomares and Pauly, 2016) in the case of invertebrates.

Ecosystem Drivers: Fishing Capacity, Nutrient Loads, and Climate Change

Data on Chioggia's fishing fleet capacity in terms of gross registered tonnage (GRT) were gathered from ISTAT (Italian National Institute of Statistics) for the period 1951–1989 and in terms of gross tonnage (GT) from the Community Fleet Register for the period 1990–2014. Total fishing capacity was expressed as GT for the entire time-series and used as a proxy for fishing pressure in the analyses because no long-term records of fishing effort were available. It was also not possible to include in fishing pressure any estimate of the technological creep due to the lack of data (e.g., CPUE from scientific surveys; Engelhard, 2016).

The monthly discharge of the Po river (m³/s), measured at Pontelagoscuro (Ferrara, Italy) for the period 1945–2014 was provided by the Regional Environmental Protection Agency of Emilia Romagna and used to calculate monthly mean annual discharge. Nutrient yearly discharge, in terms of annual input of nitrogen (NO₃) and phosphorus (PO₄) load (t/y) from the Po river, was obtained from Ludwig et al. (2009) for the period 1960–2000, and from EU FP7 PERSEUS project for the period 2001–2014.

The winter (December-March) station based index of the North Atlantic Oscillation index (Hurrell and National Center for Atmospheric Research Staff, 2016) was used to represent the large-scale climatic variability in the area in the period 1945–2014. The NAO is based on the difference of normalized sea level pressure between Lisbon (Portugal) and Stykkisholmur/Reykjavik (Iceland). Positive values of the NAO index are typically associated with stronger-than-average westerlies over the middle latitudes and wetter/milder weather over western Europe.

Finally, monthly data on sea surface temperature (SST) were downloaded from the International Comprehensive Ocean–Atmosphere Data Set (ICOADS). SST data at 13°E–45°N, 15°E–45°N, 13°E–43°N, and 15°E–43°N for the period 1957–2014 were used to compute the annual mean values.

Ecological Indicators

Data were also used to compute ecological indicators, meant to provide a holistic description of the system. Disaggregated landings per species (landed quantities in kg per year) were integrated into metrics by applying a suite of ecological indicators (described in Table 1) based on TL, thermal preference, habitat and species' size.

Table 1. Ecological indicators calculated for the time series of landings.

Metrics used included: the mean Trophic Level of landings, an indicator of food webs structure (mTL; Pauly et al., 1998); the Mean Temperature of the Catch (MTC; Cheung et al., 2013), which is the average inferred temperature preference of the species weighted by their annual catch used for evaluating the effect of sea warming on fish communities; the ratio of small pelagic fish to demersal and benthic landings (P/D ratio; de Leiva Moreno et al., 2000); the Large Species Indicator (LSI), i.e., the biomass proportion in landings of large-sized fish species (Shephard et al., 2012); the ratio between elasmobranch and bony fish in the landings (E/B ratio; Piet and Pranovi, 2005); the Primary Production Required to sustain fishery catches (PPR; Pauly and Christensen, 1995) that is the amount of energy exported from the system by landings and can be seen as the ecological footprint of fishing activities (Swartz et al., 2010); the Q-90 statistic (Ainsworth and Pitcher, 2006) that is a variant on Kempton's Q index to evaluate biodiversity (Kempton and Taylor, 1976). Landings per Unit of Fishing Capacity (LPUC) were also computed considering the yearly total gross tonnage of the fleet.

Analysis of Data

Multivariate Analysis

Landings data were aggregated into 12 functional groups, according to their taxonomical and ecological features, to reduce missing values and restrict the number of variables in order to better elucidate main changes. Invertebrates were subdivided into three groups, i.e., bivalves, cephalopods, and crustaceans (no reliable data for gastropods were available in the database and thus this group was omitted). Pelagic bony fish species were grouped into three size classes according to fish species L_max (small: L_max < 30 cm; medium: 30 cm ≤ L_max < 90 cm; large: L_max ≥ 90 cm). Analogously, three groups were used for aggregating demersal bony fish species by size. Flatfishes, sharks and “skates and rays” represented three further groups.

Functional groups' proportion in landings (%) was computed yearly in order to analyse landings composition and remove the effect of landings abundances. A fourth root transformation was applied to data to downweight the influence of predominant groups (Kaiser et al., 2000) and the similarity between every pair of years was computed using the Bray-Curtis similarity coefficient. Chronological clustering using the un-weighted pair-group average (Legendre and Legendre, 1998) was applied to identify periods with similar landings composition. Functional groups' proportion in landings was compared among periods by means of the non-parametric Kruskal-Wallis H-test (α = 0.05). Ecological indicators were compared between consecutive periods identified through cluster analysis by means of the non-parametric Mann-Whitney U-test (α = 0.05). Analysis was done with PAST v. 3.14 (PAleontological STatistics; Hammer et al., 2001).

Ecological indicators and drivers relationship was analyzed for the period 1960–2014 (for which all drivers were available) through the statistical procedure BIO-ENV from PRIMER v. 6.1.5 (global BEST-test; Clarke et al., 2008). The method consists in the computation of the correlation coefficients between similarity matrices of ecological indicators and drivers, and identifies the combination of drivers that maximizes the correlation. Indicators and drivers were normalized prior to the construction of the Euclidean distance matrices, since they represented different units of measure. The skewness and the individual correlations between drivers were explored by constructing a draftsman plot and examining the resulting Spearman rank correlations to eventually reduce redundancy and dimensionality of the data. We included all drivers in the analysis after testing the absence of highly correlated drivers (ρ ≥ 0.95).

Time-Series Analysis

Trends in total landings (Y), ecological indicators and drivers time-series were analyzed using the Mann–Kendall test and the non-parametric Sen's method was used to quantify the magnitude (slope) of the trend, using the Microsoft Excel template MAKESENS (Mann-Kendall test for trend and Sen's slope estimates) developed by Salmi et al. (2002).

Drivers time-series were further explored using the sequential t-test analysis of regime shift (STARS v. 3.4) first developed by Rodionov (2004) in order to detect potential abrupt changes in the mean. The algorithm was applied to the filtered time series calculated by removing red noise through a “pre-whitening” procedure based on the IP4 method (Rodionov, 2006) to take into account the effect of serial correlation on shift detection. The cut-off length was set at 10 years, the significance level at 0.05 and Huber's weight parameter at 1.

Results

Multivariate Analysis

Chronological clustering allowed distinguishing seven periods with different landings composition at a similarity threshold of 0.94: 1945, 1946–1954, 1955–1961, 1962–1985, 1986–1993, 1994–2008, and 2009–2014 (Figure 2). Year 1945, which represented an autonomous cluster at this similarity threshold, was associated to the first period (then becoming 1945–1954).

Figure 2. Chronological clustering (Bray-Curtis similarity, unweighted pair-group average) of the composition of landings for the years 1945–2014.

The difference among periods of functional groups' proportion in landings was statistically significant for all functional groups. Landings were dominated by small pelagics in all periods, with a percentage contribution ranging between 37% (SD: 8%; 1945–1954) and 68% (SD: 3%; 2009–2014). Cephalopods was the second most relevant group in the landings in all periods, with a contribution ranging between 10% (SD: 2%; 1955–1961) and 20% (SD: 5%; 1986–1993). In the first two periods, the third group in terms of biomass was medium pelagics that represented more than 10% of the landings [13% (SD: 4%) and 10% (SD: 6%), respectively], but then medium pelagics declined in the following periods down to 1.1% (SD: 0.4%) in recent years (2009–2014). Medium demersals ranged between 12% (SD: 2%; 1945–1954) and 6% (SD: 1%; 1955–1961). The other groups contributed with <10% in all periods. It is worth noting that skates and rays showed a sharp decline, from 3% (SD: 1%) in 1945–1954 to 0.04% (SD: 0.01%) in 2009–2014 (Figure 3). Another interesting dynamic was showed by bivalves that from 2% (SD: 1%) of the landings biomass in the first period increased to 8% (SD: 6%) in the period 1986–1993, and then sharply declined down to 1.9% (SD: 0.4%) in the most recent period (Figure 3).

Figure 3. Percentage composition of the fish community in the periods identified through chronological clustering. The confidence interval (95%) is reported in the vertical bar.

Mann-Whitney U-test results are reported in Table 2. Y, MTC, LPUC, PPR, and P/D ratio significantly increased between the first two periods, while mTL, LSI and E/B ratio significantly decreased (Figures 4, 5). Between 1955–1961 and 1962–1985 only MTC significantly increased, while mTL, LSI, and P/D ratio significantly decreased (Figure 5). In the subsequent period, mTL and LSI significantly increased, while Y, LPUC, PPR, Q-90, and P/D ratio significantly decreased (Figures 4, 5). Between 1986–1993 and 1994–2008 a significant increase of MTC, Q-90, and P/D ratio was observed, while E/B ratio further decreased (Figure 5). Finally, by comparing the last two periods (1994–2008 and 2009–2014) a significant decrease of MTC, PPR, LSI, Q-90, and E/B ratio occurred, while P/D ratio significantly increased (Figure 5).

Table 2. Significant differences (Mann-Whitney U-test, p < 0.05) of total landings (Y), landings per unit of capacity (LPUC), and ecological indicators between successive periods.

Figure 4. (A) Total landings (Y) and (B) landings per unit of fishing capacity (LPUC). Black dots, data; continuous red line, median value in each period identified through the cluster analysis; dotted red lines, interquartile range in each period.

Figure 5. Ecological indicators calculated on the basis of landings disaggregated data. Black dots, data; continuous red line, median value in each period identified through the cluster analysis; dotted red lines, interquartile range in each period. (A) mTL, mean trophic level; (B) MTC, mean temperature of the catch; (C) PPR, primary production required; (D) LSI, large species indicator; (E) E/B ratio, ratio between elasmobranch and bony fish in the landings; (F) Q-90, a variant on Kempton's Q index; (G) P/D ratio, ratio of small pelagic fish to demersal fish plus benthic landings.

The global BEST-test showed a moderate (ρ = 0.26–0.34) significant (p < 0.01) link between ecological indicators and drivers. The similarity matrices obtained with phosphorous discharge and fishing capacity correlated highest (ρ = 0.34) with the similarity matrix of ecological indicators.

Time-Series Analysis

Landings (Y) significantly increased between 1945 and 2014 (Z = 3.36, p < 0.001, slope = 58.9 t year⁻¹, Figure 4A). Conversely, landings rescaled over fishing capacity (LPUC) significantly declined in the whole period of study (Z = −3.95, p < 0.001, slope = −8.13 t GT⁻¹ year⁻¹, Figure 4B). The mTL significantly declined along the whole time series (Z = −3.38, p < 0.001, slope = −0.002 year⁻¹, Figure 5A). MTC significantly increased from 1945 to 2014 (Z = 7.27, p < 0.001, slope = 0.048°C year⁻¹, Figure 5B). PPR showed no significant trend between 1945 and 2014 (Figure 5C). LSI and E/B ratio significantly declined from 1945 to 2014 (Z = −3.21, p < 0.01, slope = −0.002 year⁻¹; Z = −8.67, p < 0.001, slope = −0.0003 year⁻¹, Figures 5D,E). The same temporal trend resulted for Q-90 (Z = −4.92, p < 0.001, slope = −0.021 year⁻¹, Figure 5F). The P/D ratio significantly increased from 1945 to 2014 (Z = 2.7, p < 0.01, slope = −0.02 year⁻¹, Figure 5G).

As regards drivers, Po river mean annual discharge did not show a significant trend. However, the time-series presented a significant negative drop in mean in 2003 and a positive one in 2008 (Figure 6A). Nitrogen load significantly increased between 1945 and 2014 (Z = 4.17, p < 0.001, slope = 1.37 t year⁻¹). The time-series presented a significant positive drop in mean in 1972 (Figure 6B). Conversely, phosphate load did not show any significant trend, but an increase in the mean value in 1972 and a significant decrease in 1988 (Figure 6C). Neither the NAO showed any significant trend, but two negative shifts in 1962 and 2009 and a positive shift between the two in 1972 (Figure 6D). SST significantly increased in the period of study (Z = 3.56, p < 0.001, slope = −0.01°C year⁻¹) and showed two positive shifts in mean in 1998 and 2011 (Figure 6E). Finally, also fishing capacity significantly increased (Z = 9.46, p < 0.001, slope = 89.88 GT year⁻¹) without any significant shift in the mean (Figure 6F).

Figure 6. Time-series of drivers of ecosystem change. Black diamonds, data; red line, regime shift. (A) Po discharge, mean annual Po river freshwater discharge; (B) NO₃, annual nitrogen flux from the Po river; (C) PO₄, annual phosporus flux from the Po river; (D) NAO, North Atlantic Oscillation index; (E) SST, sea surface temperature; (F) GT, total gross tonnage of vessels of the Chioggia's fishing fleet.

Discussion

Six periods with significantly different community structure were identified through cluster analysis of the landings composition by functional group. The most abundant group was small pelagics (mainly European anchovy Engraulis encrasicolus and European pilchard Sardina pilchardus) in all periods, followed by cephalopods in almost all periods. Major changes in community composition between periods include skate and rays and medium pelagics decline, and bivalves dome-shaped trajectory.

The increase of bivalves' proportion in the landings from 2% in 1955–1961, to 5 and 8% in the following periods, is probably linked to the introduction in the mid-1950s of the “rapido” trawl, and successively in the 1970s of the hydraulic dredge, that substantially improved fishing efficiency. However, the proportion of this functional group in the landings decreased to 3% in 1994–2008, and further decreased to 2% in 2009–2014. This sharp decrease is related to the collapse of the Mediterranean scallop (Pecten jacobaeus) in the late 1990s (data not shown) due to overfishing. In fact, between the 1960s and the 1990s, the species showed large fluctuations determined both by the intensity of fishing effort and by mass mortalities due to hypoxic conditions that occurred episodically over wide areas of the Northern Adriatic Sea (Hall-Spencer et al., 1999). A previous study (Maurizio and Castagnolo, 1986) showed that after a 6-month period of fishing, commercially sized P. jacobaeus catches collapsed by 83%, indicating the high vulnerability of the species to unsustainable exploitation, and the need for the institution of selected areas closed to trawling (Hall-Spencer et al., 1999).

Moreover, a suite of ecological indicators applied to time-series of landings provided the basis for gaining useful insights on possible causes of long-term changes. The use of multiple metrics allowed to obtain a picture of how the ecosystem has changed over the past 70 years, since different metrics emphasize distinct aspects of the underlying communities, and consequently allow disentangling the role of different drivers.

Early Changes: The Decline in Size of Fished Communities

LSI, E/B ratio, and mTL showed a significant decrease between 1945 and 2014. Among the range of ecological indicators analyzed in this paper, these are directly (LSI) and indirectly (E/B ratio and mTL) related to size. Since fishing is usually size-selective, both within and among species (Jennings, 2005), these indicators are considered to be sensitive to fishing disturbance and are expected to decrease under unsustainable exploitation.

It is worth noting that LSI does not take into account the actual size-distribution of exploited populations (that requires the availability of information of the size of caught individuals), but it is based on the L_max of species (see Table 1). Thus, the LSI changes are attributable to changes in composition of populations rather than species' size composition. Therefore, the LSI differs substantially from the Large Fish Indicator (LFI, Greenstreet et al., 2011), but it may be useful for management purposes since it requires less detailed data than LFI for highlighting effects of fishing, and can be applied also to long-term landings data series where the size structure of catches is most often not known.

The negative trend for these indicators highlights a shift toward smaller species in the community, or anyway a decrease of proportion of larger fish. Larger fish (which usually have high trophic levels) are more vulnerable to fishing and have less capacity to sustain great rates of mortality (Dulvy and Reynolds, 2002; Jennings, 2005; Myers and Worm, 2005). This phenomenon could be further exacerbated by the possible increase of small species due to the “predation release” as their predators are depleted (Bruno and O'Connor, 2005; Jennings, 2005; Myers et al., 2007).

Therefore, fishing seems to have had significant detrimental effects on the fish community composition since the first years of the time-series. Indeed in mid-1950s, the first rearrangement of the fish community structure occurred, as emerged from the cluster analysis, with a shift from large-sized high trophic level species to small planktivorous ones. Consequently, a significant negative decrease of mTL, LSI, and E/B ratio was observed.

The decline of elasmobranchs and large-sized species started even before the 1950s in the area (i.e., in the nineteenth century; Fortibuoni et al., 2010; Raicevich and Fortibuoni, 2013), but was exacerbated afterwards probably as a consequence of the industrialization of fishery and the introduction of new highly impacting fishing gears (Ferretti et al., 2013; Barausse et al., 2014). Especially linked to the decrease of skates and rays (Barausse et al., 2014; Engelhard et al., 2015), the E/B ratio decline indicates the important role of elasmobranchs as sensitive key species for detecting early signals of fisheries disturbances (Baum et al., 2003).

However, the early collapse of large-sized species is not only a consequence of the decline of elasmobranchs. In the mid-1950s medium pelagic species reached a maximum in landings (1449 kg in 1956), and dramatically declined afterwards. Medium pelagics (mainly the Atlantic mackerel Scomber scombrus) represent an important group of species for the local fish market, whose decline is hardly a consequence of changes in market demand (Meneghesso et al., 2013). Thus, it can be assumed quite conservatively that the decline in catches for medium pelagics reflects a dramatic decline of populations at sea.

Overall, the decline in mTL, E/B ratio and LSI clearly highlight that there has been a long-term fishing-down food web phenomenon (sensu Pauly et al., 1998).

Signals of Structural Changes

The Q-90 index significantly declined between 1945 and 2014, with a significant negative median-shift in the mid-1980s, indicating a reduction in the biodiversity of landings with increasing fishing impact. The change in biodiversity of landings is concurrent to the significant decrease of total landings and LPUC, suggesting that main modification of marine communities affected local fisheries production.

In the mid-1980s, regime shifts are reported for different ecosystem components in many areas, the Mediterranean Sea, the North Sea, the Baltic Sea, and the Black Sea, possibly suggesting regional effects of a larger scale northern hemispheric pattern (Conversi et al., 2010). Barausse et al. (2011) reported that a similar shift occurred also in the Northern Adriatic fish community, including also medium-high trophic level species belonging to both the demersal and pelagic habitats. Our analysis, however, while confirming the existence of the shift, do not support the hypothesis that climatic drivers played a major effect on it, and rather point to the impact of local pressures, i.e., overexploitation and nutrient loads.

Such structural changes may have altered food-web structure and impaired its functioning and resilience. Coll et al. (2008) compared North-Central Adriatic Sea food-web structure between the 1970s and 1990s, and reported a high food-web degradation regarding overexploitation of higher trophic levels and a simplification of food-web structure (lower omnivory and higher generality). Thus, it is not surprising that LPUC showed a downwards significant shift in 1986, reaching its minimum value in 2002. Only in 2014, the index recovered to values comparable to the very beginning of the time-series, when fishing technologies were markedly less developed and efficient.

Changes in the Trophic State

The significant positive trend over time in the P/D index found in the present study may depend both from eutrophication and overexploitation of resources (Libralato et al., 2004). Indeed, eutrophication and overfishing may have similar and synergistic effects on fish communities, i.e., a decline in diversity, an initial increase in productivity of benthic/demersal and pelagic food webs, then the progressive dominance of the production system by short-lived, especially pelagic species (Caddy, 1993). However, the P/D ratio in the present case seems to point out the importance of eutrophication driven events.

Pelagic fishes are generally influenced by nutrient enrichment when it stimulates the plankton production (Caddy, 1993), while demersal fishes are influenced by the dynamics of benthic community, which generally responds negatively to the conditions of excessive enrichment as shown by de Leiva Moreno et al. (2000). These authors found a mean value of P/D equal to 3.76 in the Adriatic for the historical series 1978–1988. In the present study, the P/D index had a wide dynamic trajectory, ranging between 0.31 and 4.12 with a mean value of 1.81 (± 0.84), thus confirming the Adriatic as a mesotrophic ecosystem (de Leiva Moreno et al., 2000). The index reached some peaks that can be partly related to severe anoxic/hypoxic events. For instance, anoxias in bottom waters linked to eutrophication occurred in the periods 1955–1956 and 1972–1982 (Sangiorgi and Donders, 2004), when the P/D index reached a value of ~2. Finally, the P/D index showed an increasing trend in recent years (significant upward median-shift in 2009), probably linked to a partial recovery of small pelagic species, mainly European anchovy (Carpi et al., 2015).

PPR showed a significant positive median-shift in the mid-1950s. In the mid-1980s, PPR collapsed, mainly driven by the strong reduction of small pelagic landings. As regards anchovy, its population started declining since the late-1970s, and the reasons of this abrupt reduction have been previously ascribed to climate forcing (Santojanni et al., 2006), modified inflow of Mediterranean waters in the Adriatic Sea and associated salinity changes (Grbec et al., 2002), over-fishing, increased predation of eggs and larvae by the jellyfish Pelagia noctiluca and the presence of mucilage events (Regner, 1996). PPR further significantly decreased in the last period 2009–2014.

It is worth noting that during the 1970s an increase of eutrophication occurred in the basin, lasting until the mid-1980s (Giani et al., 2012). Although nitrogen loads increased up to the mid-1970s and then maintained approximately the same values up to now (Figure 6B), phosphate peaked in the mid-1980s (Figure 6C), and its following marked decline may have resulted in limiting production. In 1985, to reduce the negative effects of cultural eutrophication (e.g., anoxia events), the nutrient load delivered to the Adriatic Sea was reduced, mainly by changing the chemical composition of soap powders (banning phosphorous from the mixture) and by improving the treatment of urban sewage and farm litter products (Decree-law No. 667 of November 25th 1985—Urgent measures to limit the eutrophication of waters).

A significant decrease of phytoplankton abundance was also observed after the 1980s, along with changes in its species composition with a shift toward smaller organisms. Such changes modified also zooplankton community (Kamburska and Fonda-Umani, 2009). This trend was a consequence of a reduction of phosphorous load, being the Northern Adriatic waters phosphorous-limited (Solidoro et al., 2009), and a decline of atmospheric precipitation and the runoff in the basin (Giani et al., 2012).

Signals of Climatic Changes

MTC significantly increased between 1945 and 2014 at a decadal rate of 0.5°C, with three significant positive median-shifts in the mid-1950s, in the early-1960s, and in the mid-1990s. A significant negative median-shift was instead observed in the last period (2009–2014), even if the median value (15.8°C) was still much higher than at the one observed at the beginning of the time-series (1945–1954: 13.2°C). Thus, overall an increasing dominance in catches of warm affinity species occurred in the landings of Chioggia, coherently with the pattern observed in other Mediterranean areas (Cheung et al., 2013; Tsikliras and Stergiou, 2014; Fortibuoni et al., 2015; Tsikliras et al., 2015). Also SST significantly increased between 1957 and 2014, with a positive shift in mean in 1998. However, SST resulted not to be one of the main ecological drivers in driving community changes in the Northern Adriatic Sea, and results from this study do not allow establishing a relationship between MTC, SST in the area and NAO.

Integrating in a Coherent Framework the Community Changes

Long-term changes in the Northern Adriatic fish community resulted to be mainly related to the impacts of fisheries and nutrient dynamics, while climate had a secondary role up to now. A coherent dynamic is obtained by scaling the PPR to N:P ratio (NO₃/PO₄) to account for the changes in the limiting factor (Solidoro et al., 2009), and comparing it with the fishing capacity through time (Figure 7).

Figure 7. Dynamics of fisheries production under changes of fishing capacity and nutrient limitation from 1960 to 2014. PPR is scaled to N:P ratio to take into account changes in the nutrient limiting factor. Solid lines represent subjective indication of average main trajectories identifiable as two alternative states of the system (state 1: line A–B; state 2: C–D–E). The mid-1980s shift is coherent with a catastrophic shift between two alternative stable states (dashed arrow from B to C) and the presence of hysteretic behavior. For going back from the actual state (state 2, lower solid line) to a state prior to the catastrophic event of mid-1980s (state 1, upper solid line) it is necessary to considerably reduce fishing capacity (left dashed line E–F, placed subjectively).

The PPR to N:P ratio increased up to the mid-1980s, showing that the exploitation development was supported by nutrient enrichment. In the mid-1980s there was an abrupt collapse after which, even for large changes in fishing capacity, PPR to N:P ratio attained much lower values (Figure 7). Notably, the mid-1980s shift is coherent with a shift between two alternative stable states, as described in Scheffer and Carpenter (2003): from the first state (line A–B in Figure 7), to a second state (line C–D–E in Figure 7). Therefore, in the mid-1980s a shift between two alternative stable states of the system likely occurred in the Northern Adriatic Sea, which is reflected also in a change in the community structure (see chronological clustering; Figure 2). The shift to a new community composition and the lower trophic potential (see previous sections) thus resulted in a much lower productive capacity (here represented by PPR to N:P ratio), even for the high fishing effort exerted in the last decades of the time series (Figure 6F). Importantly, to induce a switch back to the original state, it is not sufficient to restore the conditions present before the collapse (Scheffer and Carpenter, 2003): Figure 7 shows the hysteresis (B–C–E–F) typical to alternative stable points of complex systems. From the actual situation (state 2), to reach the same PPR to N:P ratio observed before its collapse (state 1), a relevant reduction of fishing capacity is necessary (see the hypothetical line E–F).

Although from our data it is not possible to estimate the reduction in fishing capacity necessary to switch back from state 2 to state 1, Figure 7 highlights that a reduction in fishing capacity to values lower than the ones that characterized the 1980s would be necessary, together with an intervention on nutrient loads. This result points to the fact that the current critical situation of fisheries is not solely the result of mismanagement of the fisheries sector, but is due also to other environmental policies, in particular those relative to water quality. Thus, actual management needs to account for the fact that the ecosystem is now in a different state (state 2), and the reversal of the shift of the mid-1980s implies changes in fishing capacity and nutrient balance that go far beyond those of the period of the shift.

Overcoming Limitations in the Use of Landings for Ecological Analyses

This work gives support to the possibility of using landings statistics for inferring changes in marine ecosystems through the analysis of landings composition and the application of a set of ecological indicators. Despite the intrinsic limitations of fishery-dependent data (Essington et al., 2006; Hilborn, 2007), the analysis reported here highlights that opportune indicators can produce interesting and useful assessments from landings, especially if combined with local knowledge on fisheries changes. In particular, the capability to understand relationship between landings composition to community at sea is more robust when local disaggregated landings result from multi-target and multi-gear fisheries (i.e., several distinct métiers), and when changes in fishing activities (e.g., introduction of new technologies, shift from one fishing gear to another, number and tonnage of boats) are traceable.

Being aware of the poor capabilities of landings absolute quantities to assess single species abundances (Pauly et al., 2013), we used landings composition (percentage proportion of functional groups in total landings) to detect community changes. Moreover, the use of data lumped into functional groups allowed overcoming problems of changes in aggregation detail in landings statistics.

The potential is great given that landings data are the most widespread information that can be used to analyse marine ecosystem changes in the past, and probably the cheapest information that can be collected by surveying archives and statistical bulletins and that can be used also in data poor conditions.

Conclusions

The analysis of landings from the fish market of Chioggia allowed inferring information on long-term changes in the Northern Adriatic ecosystem. Most vulnerable species (i.e., elasmobranchs and large-sized species) considerably declined at the beginning of the industrialization of fishery and continued to decline in the following decades, probably because the exploitation rates were not sustainable. Indeed, fishing capacity increased enormously during the 1960s and 1970s, i.e., larger boats, higher tonnage and engine horsepower, improved fishing gears, use of high-technology equipment (Fortibuoni, 2010). Until the mid-1980s total landings continuously increased, possibly as a consequence of the modernization of fishing fleets and of the growing cultural eutrophication. However, while landings were still increasing, the LPUC was already declining.

Later, the nutrient load (in terms of phosphorus) delivered to the Adriatic Sea decreased, thus leading to a combination of high exploitation and reduced productivity, which may well-explain the collapse in landings in the following years. Indeed, from our analysis resulted that phosphorous load (a proxy of primary production) and fishing capacity (a proxy of fishing effort) were the main drivers of change among the considered explanatory variables. Conversely, climate relate variables had a smaller impact.

It is likely that long-term effects of fishing drove significant changes in fish community structure in the Northern Adriatic Sea that were partially masked or balanced by an increase in productivity in the period of high nutrient discharge. Once productivity declined, the food-web structure was already modified and probably the resilience of the system was unpaired. The ecosystem was in a “fishing status” (sensu Jennings and Kaiser, 1998), thus reducing its recoverability from environmental driven imbalance. As regards the role played by fishing, it is expected that over time the enhanced skipper skills, adoption of auxiliary equipment and more efficient gear and materials, replacement of old vessels by new ones and upgraded engines (Damalas et al., 2014; Engelhard, 2016), resulted in an increased catching efficiency. However, since no reliable information was available for the area we did not consider the effect of the technological creep, keeping our results rather conservative with respect to fishing impacts.

This study shows that the Northern Adriatic ecosystem and ecological drivers has dramatically changed in the last decades, and thus greater knowledge of past states is crucial to set appropriate baselines for current management. Indeed, the actual crisis faced by the Northern Adriatic Sea fishery sector may be ascribed both to long-term over-exploitation and changes in nutrient load. This evidence should be considered in fishery management, for instance by rescaling the fishing capacity according to the present status of environmental parameters (e.g., trophic conditions related to nutrient discharges). Since fisheries management in Italy (as in the whole Mediterranean) is predominantly capacity/effort-based, accounting for changes in these parameters is decisive. However, because of technological creep, measuring nominal capacity and effort in conventional terms (e.g., GT, KW, days at sea) may produce estimates far from the effective fishing mortality exerted by the fleet (Damalas et al., 2014). All this may partly explain why in Italian waters the positive impact on resources expected from fishing capacity/effort reduction was lower than expected (Cataudella and Spagnolo, 2011). Thus, research on fishing power change and fine-scale analysis of fishing effort through time are highly recommended in order to understand real changes in the capacity of fishing fleets, and their potential to exploit fish stocks (Engelhard, 2016). Moreover, the results might provide evidence on the need for considering broadly the ecosystems impacts of human interventions and management actions, since actual critical state of the fisheries have been exacerbated by regulations on water quality.

Author Contributions

SL, SR, and TF conceived the study. SL, CS, and TF defined the hypotheses and methodologies. TF, SL, and CS analyzed and discussed the results with input from SR, FP, and OG. All authors were involved in reviewing and editing the manuscript, and approved it for publication.

Funding

The authors acknowledge partial support from EU FP 7 project PERSEUS (Policy-oriented Marine Environmental Research for the Southern European Seas, http://www.perseus-net.eu/) and the Italian RITMARE Flagship Project—The Italian Research for the Sea (http://www.ritmare.it) coordinated by the Italian National Research Council and funded by the Italian Ministry of Education, University and Research within the National Research Program 2011–2013.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge Paolo Lazzari (National Institute of Oceanography and Experimental Geophysics, Italy), Virna Meccia (Max Planck Institute for Meteorology, Germany), Cristina Mazziotti and Giovanna Pirretti (Regional Environmental Protection Agency, Italy) for helping us with the collection of ecosystem drivers data (nutrients load, climatic data, Po river discharge), and Carlotta Mazzoldi (University of Padua, Italy) for providing us landings data. This study was benefited from discussions with colleagues of the ICES Working Group on the History of Fish and Fisheries. This article is based upon work from COST Action IS1403, Oceans Past Platform (OPP), supported by COST (European Cooperation in Science and Technology).

References