Catch data (abundance and biomass) of five Macrouridae species (C. caelorhincus, H. italicus, N. aequalis, N. sclerorhynchus and T. scabrus) were analysed from a total of 22 annual cruises of the MEDITS survey programme (1994-2015). The study area covered the northern Mediterranean, from the northern Alboran Sea to the Aegean Sea, and comprised six countries (Spain, France, Italy, Albania, Montenegro and Greece) and 15 different geographical sub-areas (GSAs) (Fig. 1). Although there are gaps in the time series of some GSAs (2: 1995-1999, 2001-2005, 2010; 5: 1994-2006; 8: 2002; 20, 22 and 23: 2002, 2007, 2009-2013, 2015) because MEDITS surveys were not conducted in several years for specific reasons, they have been included to cover the largest geographical range. However, the northern Adriatic Sea (GSA 17) was excluded from this study due to the absence of Macrouridae species in its shallow seabed features.

For each haul, the number and weight of individuals belonging to the Macrouridae species were standardized to one-hour towing to calculate both species abundance (number of individuals per hour of towing time [ind h⁻¹] and biomass [kg h⁻¹]). The mean fish weight was obtained for each species by dividing biomass by abundance. Data were pooled in a matrix of species abundance, biomass and mean weight according to each haul (species vs haul). Frequency of occurrence of each species (F) was expressed as a percentage and was calculated, for the whole area, as the ratio of the number of occurrences of a species to the total number of hauls and, for each GSA, as the ratio of the number of occurrences of a species to the total number of hauls per GSA. Since the abundance of the macrourids studied was negligible at depths shallower than 200 m, further analysis was confined to bathymetric strata deeper than 200 m.

The centre of gravity (COG) was computed for abundance data in order to describe the mean bathymetric distribution of each species and to indicate the depth interval in which the species reaches its maximum abundance (Stefanescu 1991). For this analysis, the bathymetric range (200-800 m) was divided into six 100-m strata, and COG was calculated as COG=(x₁+2x₂+3x₃+…+nx_n)/∑x_i, where x_i represents the calculated mean abundance value of species x in stratum i.

Generalized additive models (GAMs) were applied to assess the bathymetric, geographic and temporal effects on abundance, biomass and mean weight of species by haul. Year was considered as a factor. A one-dimensional smoother was used to investigate the bathymetric effect, while a two-dimensional smoother was used to account for the geographic effect, combining latitude and longitude. The logarithmically transformed values (log[x+1]) of abundance and biomass were used in order to ensure a Gaussian distribution of the residuals. For the selection of the best model for each response variable, minimization of the Akaike information criterion (AIC) was applied. The complete applied model for each response variable (RV, log-transformed values of abundance and biomass, and mean weight) and each Macrouridae species was as follows:

RV=factor (Year)+s(Depth)+g(Latitude, Longitude)+ɛ,

with s and g as the univariate and the bi-variate smoothers, respectively, and ɛ denoting the Gaussian error term. The package mgcv in R (http://www.r-project.org) was used to apply the GAMs (Wood 2006).

Dynamic factor analysis (DFA) was used to identify common trends in standardized abundance data series among GSAs. DFA is a multivariate analysis belonging to dimension reduction techniques; it is designed for relatives in which a set of time series is modelled as a linear combination of underlying common trends, factor loadings and error terms to explain the temporal variability. Correlation of observation errors was modelled using different error matrices (Zuur et al. 2003). For all species, the correlations of observation errors were fitted to all possible model structures in the time series, including one to three common trends (Keller et al. 2017). The corrected AIC (AICc) was used to measure for goodness of fit having the best model the lowest AIC (Zuur et al. 2003). Model implementation was done in R (version 3.3.3), using the multivariate autoregressive state-space (MARSS) package.

Scientific names for species followed the nomenclature of the World Register of Marine Species (WoRMS 2017).

Coelorinchus caelorhincus was present in every GSA sampled and was caught on 4107 of the 8289 hauls undertaken on the continental slope (F=50%). The highest F values were observed in the northern Alboran Sea and Corsica (GSAs 1 and 8), followed by Alboran Island and the Gulf of Lions (GSAs 2 and 7) (Table 1).

GAM analysis showed a statistically significant effect of depth and latitude-longitude on log-transformed indices of abundance, biomass and mean weight of the species. The best models, explaining 40.3%, 42.9% and 35% of the deviance (abundance, biomass and mean weight of the species, respectively), included depth, latitude-longitude and year as a factor (Table 2).

The bathymetric effect showed a non-linear pattern, with biomass and abundance peaking at 400 to 500 m depth. As depth increases, the mean weight of the species also increases, with a maximum value between 500 and 600 m depth. Abundance and biomass were highest in the northern Alboran Sea, Alboran Island, the Gulf of Lions and southern Sicily (GSAs 1, 2, 7 and 16). Maximum mean weight was registered in the Gulf of Lions, Corsica and the Ligurian Sea (GSAs 7, 8 and part of 9) (Fig. 5).

The DFA model for C. caelorhincus abundances showed a common and general rise in abundance in the whole region, with the highest recorded values corresponding to the last years of the time series. An exception to this pattern was found for the Balearic Islands, Sardinia and Crete (GSAs 5, 11 and 23) (Fig. 6).

This species was present in every GSA sampled and recorded on 4934 of the 8289 hauls undertaken on the continental slope (F=60%). The highest F values were observed in the southern and central Tyrrhenian Sea and southern Sicily (GSAs 10 and 16) and the lowest in the northern Alboran Sea and Alboran Island (GSAs 1 and 2) (Table 1).

The best GAM models explaining 51.3%, 54.2% and 7.4% of the deviance (for abundance, biomass and mean weight of the species, respectively) included depth, latitude-longitude and year as a factor (Table 2).

Abundance and biomass peaked between 400 and 600 m, whereas the mean weight showed a positive relationship with depth. An increase in abundance and biomass of the species was detected over the central Mediterranean Sea, while no clear geographical effect was observed for mean weight (Fig. 5).

For H. italicus, DFA results showed one common non-linear trend with two periods of generally higher values. However, this trend was not geographically homogeneous and had a clear geographic segregation: positive loadings for most eastern areas, except for the Aegean Sea (GSA 22), and negative (i.e. generally opposite trend) for western areas (Fig. 6).

This species was absent in Corsica, the southern Adriatic, the western and eastern Ionian Sea, Aegean Sea and Crete (GSAs 8, 18, 19, 20, 22 and 23) during the period analysed. It was caught on 842 of the 8289 hauls (F=10%). The highest F values corresponded to Spanish GSAs and were null or very low in the other areas (Table 1).

The best GAM models explaining 53%, 51.3% and 26.8% of the deviance (for abundance, biomass and mean weight of the species, respectively) included depth, latitude-longitude and year as a factor (Table 2).

Abundance, biomass and mean weight increased with depth. The abundance and biomass of N. aequalis increased in the western areas, mainly in Spanish waters and particularly in the northern Alboran Sea and Alboran Island (GSAs 1 and 2). However, the geographical pattern of mean weight was quite heterogeneous, with noteworthy small-scale differences. Low mean weight values were observed along the Spanish coasts, with the exception of the high values recorded for GSAs 1 and 2 (Fig. 5).

For this species, DFA showed the lowest values during the first years and the highest values around 2006. Positive loadings for most western areas suggest a general positive relationship of the abundances with this trend. Negative loadings were obtained for the Sardinian Sea and southern Sicily (GSAs 11 and 16) (Fig. 6).

This species was not recorded in the northern Alboran Sea, Alboran Island, the Balearic Islands, northern Spain and Crete (GSAs 1, 2, 5, 6 and 23). It was caught on 2927 of the 8289 hauls (F=35%), and the highest F values were recorded in the western Ionian Sea (GSA 19), followed by the southern Sicily and the southern Adriatic Sea (GSAs 16 and 18) (Table 1).

The best GAM model explaining 54.5%, 55.7% and 17.9% of the deviance (for abundance, biomass and mean weight, respectively) included depth, latitude-longitude and year as a factor (Table 2).

Abundance and biomass displayed sharp increases from 500 m depth, with maximum values around 700 m and a further decrease to 800 m. The species mean weight increased slightly from 300 to 700 m and sharply below this depth range. The species increased its abundance in the central and eastern Mediterranean areas, but showed an opposite pattern when the mean weight was analysed (Fig. 5).

The DFA showed one non-linear but general increasing trend, with relatively stable values from 1998 to 2006. All areas showed positive factor loadings, with the exception of Corsica (GSA 8) (Fig. 6).

The species was not recorded in the southern and central Tyrrhenian Sea, southern Sicily, the eastern Ionian Sea, the Aegean Sea and Crete (GSAs 10, 16, 20, 22 and 23). It was recorded on 1035 of the 8289 hauls (F=13%). The highest F values were reached in Alboran Island, the Gulf of Lions and the northern Alboran Sea (GSAs 2, 7 and 1) (Table 1).

The best model explaining 39.3%, 41.2% and 48.5% of the deviance (for abundance, biomass and mean weight, respectively) included latitude-longitude, depth and year as a factor (Table 2).

Species abundance and biomass rose from 300 to 700 m depth, with a noteworthy increase from 700 to 800 m depth. The mean weight remained stable to 500 m and increased continuously beyond this depth to 800 m. This species was mainly present in the western Mediterranean, particularly in the northern Alboran Sea, Alboran Island and the Gulf of Lions (GSAs 1, 2 and 7), with a consistent pattern of highest mean weight of the species in these areas (Fig. 5).

Overall, increasing trends were detected by DFA models, the highest values being reached in 2006. After this year, the lowest values were reached in 2008, but they increased again thereafter. All areas showed positive factor loadings, with the exception of the southern Adriatic Sea (GSA 18), with the highest values recorded in the westernmost Mediterranean areas (GSAs 1, 2 and 5) (Fig. 6).

Depth and geographical distribution trends in abundance, biomass or mean weight differed among the five macrourid species. Abundance and biomass of C. caelorhincus increased with depth and peaked at 400 to 500 m. A similar trend was seen for H. italicus, but the peak was reached at 400 to 600 m. These results concur with those of studies carried out in the Catalan Sea (Massutí et al. 1995), in the western Ionian Sea (D’Onghia et al. 2000) and along the Iberian Mediterranean margin (Moranta et al. 2007). Massutí et al. (1995) reported a progressive decrease in H. italicus below 600 m depth. However, this species has been collected at a maximum depth of 1202 m in the central-western Mediterranean (South Sardinian) (Follesa et al. 2011). Regarding N. sclerorhynchus, abundance and biomass increased with depth, peaking at 500 to 700 m. Previous studies had shown this species reaching depths of 1500 m in the Ionian Sea (D’Onghia et al. 2004a, b) and of 1598 m in southern Sardinia (Follesa et al. 2011). Abundance and biomass of N. aequalis increased with depth, while for T. scabrous they rose abruptly below 700 m. The maximum depths reached by N. aequalis and T. scabrus were 1500 and 1200 m, respectively in the Balearic Sea (D’Onghia et al. 2004a, b), 1500 m for both in the Catalan Sea (Tecchio et al. 2013), and 1598 and 1202 m, respectively, in the southern Sardinian Sea (Follesa et al. 2011). In the Ionian Sea, T. scabrus was caught as deep as 1000 m (D’Onghia et al. 2004a, b).

In marine environments, depth is one of the main factors influencing biological distributions and is generally coupled with the combination of the gradients that co-occur with depth, affecting the biology and physiology of marine organisms and the ecological interactions between taxa (e.g. Rex 1977, Carney 2005, Drazen and Haedrich 2012). On the Mediterranean continental slope, given the prevailing conditions of relative stability below 200 m (Hopkins 1985), several authors suggest that depth-related changes are due to biological causes (e.g. Stefanescu 1991, Cartes and Carrassón 2004), such as the dispersal capabilities of early developmental stages and trophic causes (Rex 1977, Grassle et al. 1979). The existence of different niche dimensions enables both coexistence and habitat substitution of morphologically neighbouring species of Macrouridae along environmental gradients (e.g. depth). All Mediterranean macrourids are generalist feeders (Macpherson 1979), which probably constitutes an adaptive advantage in low-productivity deep-water environments inhabited by macrourids (Madurell and Cartes 2006). Their ecological segregation is maintained by a combination of differential depth distributions and feeding habits (Carrassón and Matallanas 2002). Food partitioning among coexisting macrourid species is usually related to differences in their bathymetric distribution and in their morphology (Macpherson 1979, Mauchline and Gordon 1984, Carrassón and Matallanas 2002). The different head morphologies of macrourids are correlated with their various diets and foraging tactics. Furthermore, the bathymetric centres of distribution also tended to be different, reinforcing the dietary differences as a means of niche separation (Mauchline and Gordon 1984). Macpherson (1979) showed that rates of competitive exclusion were low among four species of macrourids (C. caelorhincus, H. italicus, N. aequalis and T. scabrus) on the upper-middle slope (200-800 m depth) of the western Mediterranean Sea. Even so, H. italicus overlaps considerably with the other species, probably due to its more pelagic habitats. Copepods, amphipods and other pelagic crustaceans form the main fraction of the diet of H. italicus. The diets of N. aequalis and C. caelorhincus consist largely of polychaetes, isopods, amphipods, mysids and decapod crustaceans, while T. scabrus feeds heavily on decapods (Macpherson 1979).

The five species studied showed a general pattern of mean weight increasing with depth. The bigger-deeper phenomenon has been previously described for Mediterranean macrourids as a well-defined rule (Massutí et al. 1995, D’Onghia et al. 2000, Moranta et al 2007). Macpherson (1979) detected positive size-depth correlations in macrourids, as well as an increase in the mean size of prey with the body size of these species. Indeed, T. scabrus shows prey changes, since it mainly consumes Calocaris macandreae at around 600 m depth, and its bathymetric distribution seems to be strongly linked to that of its preferred prey. However, below 1000 m, larger T. scabrus consume more suprabenthic-bathypelagic prey (Carrassón and Cartes 2002, Cartes et al. 2009).

Emig and Geistdoerfer (2004) stated that species occurring in both the western and eastern Mediterranean basins are always found at lower depths in the eastern one. In this study, the COG shows that C. caelorhincus, the macrourid species analysed with the most homogeneous distribution throughout the Mediterranean Sea, appears to be found deeper in some of the eastern areas. However, this pattern has not been identified for the other species.

Different relationships were detected between the abundance and biomass of the five macrourid species and different geographical areas of the Mediterranean basins. C. caelorhincus and H. italicus occurred in the entire area studied, unlike the other species analysed, but their distribution patterns were different. The shallowest species, C. caelorhincus, with relatively high catch levels in most areas, seems to be the best adapted to the environmental differences between the western and eastern Mediterranean basins. However, differences among areas for H. italicus were evident, since it was more abundant in central and easternmost Mediterranean areas. This is in agreement with D’Onghia et al. (2000), who reported H. italicus in 80% of trawl hauls carried out in the Ionian Sea, and with Moranta et al. (2007), who recorded low catches along the Spanish Mediterranean coast. They were attributed to the more pelagic feeding habits that might make this species less catchable by bottom trawls. Since our study was carried out with the same gear throughout the Mediterranean Sea, differences between areas may be due to other causes. This species appears to favour warmer water masses, since it is commonly found in the tropical Atlantic and western Indian Ocean (Froese and Pauly 2017). Indeed, Sobrino et al. (2012) identified H. italicus as one of the most abundant species in African Atlantic waters.

Large differences between the western and eastern Mediterranean basins were found for N. aequalis, T. scabrus and N. sclerorhynchus. Basically, N. aequalis and T. scabrus were reported in the western part of the Mediterranean, with low or no catches in other areas. By contrast, N. sclerorhynchus was more abundant in the central and eastern parts. These species also seemed to display latitudinal differences throughout the Atlantic Ocean, with a more northward distribution for N. aequalis and T. scabrus than for N. sclerorhynchus, which has a more tropical Atlantic distribution (Froese and Pauly 2017). Differences between the western and eastern Mediterranean Sea in N. aequalis and N. sclerorhynchus were previously reported by D’Onghia et al (2004a, b) in the Balearic Islands and Ionian Sea. These authors argued that the overlapping diet of the two species may be responsible for their competitive exclusion, although the presence of both species on the western side of the Mediterranean Sea has long been known (e.g. Rannou 1975, 1976). In our study, N. aequalis was more abundant in the Alboran Sea (northern Alboran Sea and Alboran Island) than in the remaining areas, with a sharp decrease in the Catalan Sea and Balearic Islands. None of these areas reported any N. sclerorhynchus. Both species were caught in the Gulf of Lions, the Ligurian, Tyrrhenian and Sardinian seas, and southern Sicily, although catches of N. aequalis were very low in these areas. Likewise, both species were reported off the southern Sardinian Sea (Follesa et al. 2011). Therefore, these areas between the western and eastern Mediterranean Sea could be transition zones for the species, where N. aequalis begins to disappear. In Atlantic waters, both species were found by Sobrino et al. (2012), while Serrano et al. (2011) did not report N. sclerorhynchus. However, D’Onghia et al. (2004b) also pointed out the misidentification of the two species as the underlying reason for different species distribution. This prospect should also be considered in this study.

According to our results, T. scabrus was also more abundant in western Mediterranean areas, showing higher values of abundance, biomass and mean weight in the Alboran Sea and Gulf of Lions than in the remaining areas. It is widely known that these two areas display the highest primary production values for the Mediterranean Sea (Vargas-Yáñez et al. 2010). This might indicate some species’ preference for more productive systems that may facilitate a more efficient bentho-pelagic coupling and energy transfer efficiency.

Overall, in this study, the high recorded values of abundance, biomass and mean weight were recurrent for some species in the Alboran Sea and Gulf of Lions. Local hydrographic features and topographic differences greatly influence the spatial variability of the environmental parameters within each sub-basin (Tselepides et al. 2004). Regarding the Alboran Sea, Moranta et al. (2007) have already found differences in abundance and body size for Macrouridae species along the Iberian Mediterranean coast, pinpointing the particular oceanographic conditions of the Alboran Sea as a possible cause. These authors finally concluded that the impact of fishing exploitation could mask the effect of abiotic factors since, unlike the other Iberian Mediterranean areas, the open slopes of the western Alboran Sea have remained almost unexploited below depths of 500 m (Gil de Sola 1993). Other deep-sea fish species show differences between Alboran and the rest of the Mediterranean. Galeus melastomus and Chimaera monstrosa are more abundant in the Alboran Sea (Massutí and Moranta 2003). Cabo de Gata (located further east of the Alboran Sea) is a boundary for Galeus atlanticus, which is absent in the rest of the Mediterranean (Rey et al. 2010). However, these species are also abundant in Atlantic waters, on the northern Iberian slope (Serrano et al. 2011). The particular geomorphological and associated oceanographic conditions of the Alboran Sea make this region a distinct transition zone between Mediterranean and Atlantic fauna (Abelló et al. 2002).

Regarding inter-annual trends, abundance of macrourid species showed no decreasing trends over the 22 years analysed, despite the fact that these species constitute a large fraction of the discards from the deep-water bottom trawl fishery (Moranta et al. 2000, D’Onghia et al. 2003). In agreement with this, Granger et al. (2015) found that the diversity of Mediterranean demersal fish species did not decrease during the period 1994-2012. Fishery pressure causes changes in the composition of marine communities over short time scales, especially in vulnerable habitats such as deep-sea grounds. Decreases in the catch per unit effort and size of exploited deep-sea fish have occurred just a few decades after the commencement of exploitation (Cartes et al. 2013). Granger et al. (2015) suggested that human pressure had probably already impacted fish diversity prior to 1994. Overall macrourids are known to have adapted to various ranges of food availability by employing generalist and opportunistic foraging strategies (Pearcy and Ambler 1974, Macpherson 1979, Mauchline and Gordon 1984). In addition, they can have a continuous spawning span in the Mediterranean Sea (Massutí et al. 1995, D’Onghia et al. 1996, 1999, Fernandez-Arcaya et al. 2013) and can be found in habitats less impacted by fishing activities, such as submarine canyons and cold-water coral communities, which can act as shelters from trawling (Tursi et al. 2004, D’Onghia et al. 2010, Fernandez-Arcaya et al. 2013). Most macrourid species find refuge in deep waters since they exceed the depth range sought by fishing in the Mediterranean (D’Onghia et al. 2004a, b, Follesa et al. 2011, Tecchio et al. 2013).

Inter-annual abundance increasing trends were detected for C. caelorhincus and N. sclerorhynchus in almost every area studied. The two deepest species, N. aequalis and T. scabrous, showed similar inter-annual trends, with abundances peaking in 2006 and fluctuating thereafter. No consistent increasing or decreasing trends, with differences between western and eastern Mediterranean areas, as reported above, were observed for H. italicus. Trophic factors are the most likely link between climatic and physical oceanographic processes and the abundance and fisheries yields of top predators (Cartes et al. 2009). Cartes et al. (2009) compared recent (2007 and 2008) and former data (from 1988 to 1992), finding a dominance of planktonic suprabenthic feeders (e.g. fish such as H. italicus) in the former time series (1988-1992) and a higher abundance of benthic feeders in 2007, coinciding with changes recorded in the North Atlantic Oscillation index.

The Mediterranean Sea is a complex ecosystem with contrasting regions in terms of productivity (Tecchio et al. 2011), seafloor topography and hydrography (Millot 2005). In an ecosystem with such pronounced regional differences, animal populations are prone to exhibit patchy distributions due to different habitat conditions (Keller et al. 2017), such as the spatial distribution throughout the Mediterranean of the five macrourid species described and discussed here for the first time. Our results indicate that both depth and latitude-longitude are determinant factors of their presence and abundance. Additionally, no consistent inter-annual trends in abundance of the species were found over the 22-year-study period. Further studies taking into account spatial and temporal environmental factors, including anthropogenic variables such as fishing activity throughout the Mediterranean Sea, are required in order to elucidate the role of the variables associated with fish species distribution and their temporal changes.