Biological and environmental samples were collected during several cruises in the framework of the ESF EuroDEEP BIOFUN project (“Biodiversity and Ecosystem Functioning in Contrasting Southern European Deep-sea Environments: from Viruses to Megafauna”), whose main aim was to explore possible links between biodiversity patterns and ecosystem functioning in relation to environmental conditions, through analyzing organisms ranging from viruses to megafauna.

Sampling was conducted at selected deep-sea sites along a gradient of increasing oligotrophy from west to east—from the Galicia Bank in the Atlantic Ocean to the Levantine Basin—at depths of 1,200, 2,000, and 3,000 m. Sediment samples were collected in 7 open-slope areas: 1 in the North Eastern Atlantic (ATL), 3 in the Western Mediterranean basin (wM1, wM2, and wM3), and 3 in the Central-Eastern Mediterranean basin (c-eM1, c-eM2, and c-eM3) (Figure 1). In each area, 3 stations were sampled at 3 different depth ranges: upper bathyal (1,200 m), mid-bathyal (1,800–1,900 m), and lower bathyal (2,400–3,000 m). Since the sample at the lower bathyal depth (c-eM1) could not be collected, this station was replaced with one at a depth of 2,120 m (Table S1).

The Atlantic sampling area was on the Galicia Bank, a seamount situated on the Iberian margin about 200 km off the Galician coast. The Bank, whose summit is found at a depth of 620 m, is separated from the shallower parts of the continental margin by the Galicia Interior Basin (Pape et al., 2013), whose approximate depth is 3,000 m. Water current velocities on top of the seamount range from 5 to 30 cm s-1 (Pape et al., 2013), which is sufficient to influence organic matter deposition. Such currents result in very low phytopigment and biopolymeric organic carbon concentrations and in the presence of coarse sediments at 1,200 m depth (Table S1). The deep Mediterranean Sea is a highly oligotrophic environment, where nutrient depletion in surface waters combines with high water temperature, promoting the degradation of sinking organic matter (Sardà et al., 2004). The Mediterranean basin is characterized by a trophic gradient between the less productive eastern basin and the more productive western basin (D'Ortenzio and Ribera d'Alcalà, 2009; Giovannelli et al., 2013; Table S1), where the greater nutrient inputs are due to river runoff, the inflow of Atlantic surface water, and the outflow of relatively nutrient-rich Levantine Intermediate Water through the Strait of Gibraltar (Bergamasco and Malanotte-Rizzoli, 2010). The Mediterranean open-slope systems selected for the present study are characterized by a high percentage of silt sediment at all depths. Bottom-water temperature and salinity increase significantly from west to east, with values ranging respectively from 13.1°C and 38.5°C in the wM area to 14.7°C and 38.8°C in the Eastern Mediterranean basin (Baldrighi et al., 2014; Table S1). At each station, sampling involved 3 independent box-corer deployments for macrofauna (entirely sieved) and 3 deployments for meiofauna, microbial components, and environmental variables (n = 3 replicate samples for each parameter). Macrofauna samples were collected with a cylindrical box-corer with an internal diameter of 32 cm at all stations except for ATL, wM2 and c-eM2, where its internal diameter was 50 cm, accounting for a surface area of 803.84 cm² and 1,962.5 cm², respectively.

Data of macrofaunal abundance, biomass and diversity as well as data of environmental parameters characterizing the Mediterranean basin from west to east have previously been shown in Baldrighi et al. (2014). Data of meiofauna abundance, structural and functional diversity from the Mediterranean basin have been presented in Baldrighi and Manini (2015) and the relationship with macrofaunal structural and functional diversity has been investigated. In the present study, the same dataset has been used, with the addition of new data from the Atlantic area (i.e., macrofaunal diversity, benthic biomass, environmental features) and ecosystem functioning and efficiency proxies. We adopted a new approach of investigation aimed at finding BEF relations. A complete species list is provided in Table S2.

For grain size, biochemical composition of organic matter, and microbial component analysis, sediment subsamples from each box-corer were collected using plexiglass cores with an internal diameter of 3.6 cm. The top 1 cm of a sub-sample from each box-corer was collected and frozen at −20°C to analyze chlorophyll-a, phaeopigment, and organic matter content. Replicates of about 1 ml wet sediment were fixed using buffered formaldehyde and stored at +4°C until processing for total prokaryotic abundance and biomass determination (Giovannelli et al., 2013). The top 20 cm was preserved at +4°C for grain size analysis. For meiofauna analysis, sediment from each box-corer was collected using a plexiglass tube with an internal diameter of 3.6 cm and immediately fixed in 4% buffered formalin and Rose Bengal. Once in the laboratory, only the top 5 cm was sieved through a 300-μm and subsequently a 20-μm mesh sieve. Meiofaunal samples were obtained from 6 of the 7 areas (meiofauna data not available from c-eM2 slope area). For macrofauna analysis, the top 20 cm of sediment from each box-corer and the overlying water were sieved through a 300-μm mesh sieve to retain all macrobenthic organisms (here considered sensu lato, including nematodes, copepods, foraminifera, etc.) (Baldrighi and Manini, 2015). The residue left in the sieve was immediately fixed in buffered formalin solution (10%), and stained with Rose Bengal.

For grain size analysis, aliquots of fresh sediment were sieved through a 63-μm mesh. The two fractions (>63 μm, sand; <63 μm, silt and clay) were dried in an oven at 60°C and weighed. Data were expressed as a percentage of the total sediment dry weight (Baldrighi et al., 2014).

Chlorophyll-a and phaeopigments were determined by standard techniques (Danovaro, 2010). The sum of chlorophyll-a and phaeopigment concentrations was designated as chloroplastic pigment equivalents (CPE). Total phytopigment concentrations were converted to carbon (C) equivalents using a conversion factor of 40 (De Jonge, 1980) and expressed as mgC g⁻¹. Biopolymeric organic carbon (BPC) was calculated as the sum of the carbon equivalents of carbohydrates, proteins, and lipids (Fabiano et al., 1995). The protein to carbohydrate (PRT: CHO) ratio was then calculated and used as a descriptor of the nutritional value of sediment organic matter. PRT: CHO ratios >1.0 indicate relatively high quality and high food availability (Pusceddu et al., 2010).

The total prokaryotic number (TPN) was determined by acridine orange staining (Luna et al., 2002) and analyzed using epifluorescence microscopy (magnification, 1,000 x). Total prokaryotic biomass (TPB) was estimated using an ocular micrometer; prokaryotic cells were assigned to different size classes (Fry, 1990) and converted to biovolumes based on the assumption that they had an average carbon content of 310 fgC μm⁻³ (Fry, 1990). For meiofauna analysis, the fraction remaining on the 20-μm sieve was resuspended and centrifuged using Ludox HS40 (Danovaro, 2010). Meiofaunal biomass was estimated as described by Danovaro (2010). Macrofaunal organisms were sorted under a stereomicroscope and identified to the lowest possible taxonomic level. For each species, the total number was calculated and the wet weight biomass measured and expressed respectively as abundance and biomass per m². Wet biomass (g wet weight m⁻²) was converted to ash-free dry weight and organic carbon content using standard conversion factors (Rowe, 1983).

The diversity and functional diversity of the macrobenthic community were analyzed in each open-slope area (Table S3). Diversity was assessed by 3 measures: richness of macrofaunal higher taxa (no. of taxa), species richness (SR) in each box-corer sample, and expected number of species in each sample of 50 individuals (ES_(n)) (Table 1). The SR index and the rarefaction index were both calculated, because they are sensitive to low abundances (Carney, 1997); in particular, the latter measure provides a useful tool to compare species richness among samples with different total abundances or from different surface areas (Soetaert and Heip, 1990). On the other hand, SR incorporates macrofaunal abundance.

Functional diversity was assessed using 4 indices (Table 1): trophic diversity (Θ⁻¹), expected number of deposit feeders (EDF₍₃₀₎), expected number of predator species (EPR₍₂₀₎), and bioturbation potential (BP) (Baldrighi and Manini, 2015; Quéiros et al., 2015). The trophic diversity index was included, because affinity to a given trophic group entails a discrete selection of food sources, hence a different ecological role (Dunne et al., 2002), and feeding mechanisms are held to be central processes structuring marine ecosystems and providing a strong link between species and ecosystem functions (Pearson and Rosenberg, 1987). The number of predator species feeding on macrobenthic organisms was calculated, because the number of species at the top of the benthic food web may reflect a higher functional diversity of the entire benthic assemblage (Ngai and Srivastava, 2006). The number of deposit feeders was included, because they were the most numerous trophic group in our macrofauna samples, and deposit feeding is one of the most effective feeding strategies in environments with generally meager food sources, such as the deep sea (Gage and Tyler, 1991). Bioturbation activity on and in the seafloor, is closely related to the way organisms feed and move (Welsh, 2003) and it is assumed to have an influence on benthic processes (Belly and Snelgrove, 2016). Indeed, we considered the macrofaunal bioturbation activity to be an important functional trait. In the present work, we estimated the macrofaunal bioturbation potential according to the formula reported in Quéiros et al. (2015):

where Bi and Ai are the biomass and abundance of species/taxon i in a sample, Mi is mobility and Ri is sediment reworking. Trait scores of different macrofaunal organisms were assigned following the table of values suggested by Queirós et al. (2013), based on an extensive review of published material, and by Shields and Hughes (2009).

Deep-sea ecosystem functioning was estimated as benthic faunal biomass (mgC m⁻²), considering total benthic biomass (the sum of prokaryotic, meiofaunal, and macrofaunal biomass) and the biomass of the functional group of macrobenthic predators (Table S4). Standing biomass is a metric of ecosystem functioning (e.g., Mora et al., 2011), and a reduction in the size of the predator population may exert effects that go beyond top-down control, thereby affecting cross-system connectivity and ecosystem stability (McCauley et al., 2015). Four indicators were used to assess ecosystem efficiency (Table 1): (i) the ratio of macrofaunal biomass to biopolymeric carbon content (MBM: BPC), a measure of the ability of the system to channel detritus to higher trophic levels (Danovaro, 2012); (ii) the ratio of macrofaunal biomass to prokaryotic biomass (MBM:TPB); (iii) the ratio of macrofaunal biomass to meiofaunal biomass (MBM:MEB); and (iv) the ratio of meiofaunal biomass to prokaryotic biomass (MEB:TPB). Several deep-sea macrobenthic organisms are considered to be deposit feeders because they ingest large amounts of sediment together with detritus, prokaryotes, and meiofauna (Gage and Tyler, 1991). It has been suggested that up to 24% of total bacterial production is grazed by macrofauna, and the meiofauna are assumed to be an important link between smaller (e.g., bacteria) and larger organisms (e.g., macrofauna) (Van Oevelen et al., 2006). The MBM:TPB and MBM:MEB ratios are thus measures of the energy transferred from lower to higher trophic levels based on the hypothesis that macrofauna prey upon microbial and meiofaunal components. The MEB:TPB ratio was included because meiofauna production greatly contributes to total secondary production, particularly in deep-sea systems, where smaller organisms take over in term of standing stock (Rex et al., 2006).

The presence and nature of BEF relationships can be affected by the spatial scale of the investigation and by the environmental factors acting on each scale (Gamfeldt et al., 2014). The presence of a BEF relationship was investigated on two scales: (i) a large spatial scale encompassing the entire dataset, where the data from all 3 research areas, 7 open slopes, and 3 depths were pooled and subjected to statistical analysis; and (ii) a basin scale, where the dataset from each sampling area (ATL, wM, and c-eM) was analyzed separately. The relationships among BEF, biodiversity and ecosystem efficiency were estimated by a linear model (in the form y = a+bx), a power model (y = a+x^b), and an exponential model (y = e^a+bx). These models are considered to be the best tools to describe BEF relationships in different deep-sea environments (Cardinale et al., 2007; Danovaro et al., 2008a; Lefcheck and Duffy, 2014). Statistical analysis was performed using R-cran software (http://www.R-project.org). Map plots were drawn using Ocean Data View (Schlitzer, 2011). Relationships between variables were tested using linear and non-linear regression. After fitting the three models to the experimental data, the distribution of the residuals, r², and the Akaike Information Criterion (Akaike, 1974) were used to discriminate the best-fitting model, as appropriate. Model fitting was performed for the large scale, (i.e., the entire dataset), and the basin scale, [i.e., each sampling area (ATL, wM, and c-eM)]. Distance-based multivariate regression analysis with forward selection (DISTLM) (Anderson, 2004) was used to account for the potential effect of environmental variables on the BEF relationship. The effects of depth, longitude, temperature, grain size, and food source quantity and quality were included as covariates in the analyses. In brief, environmental covariates were assessed independently against biodiversity and ecosystem functioning (BEF) variables using a multivariate linear regression model before testing for BEF relationships. P-values for each covariate were obtained with 4,999 permutations of residuals under the reduced model (Anderson, 2001).

In addition, the collinearity between the descriptors of taxonomic and functional diversity was assessed using the Variance Inflation Factor (VIF) test (O'Brien, 2007). Based on the rule of 10, which is commonly used in statistical analysis (O'Brien, 2007), VIF values were greater than 3.0, indicating the lack of a strong correlation between our diversity indices. However, use of a threshold value of 2.50 (which corresponds to an R² > 0.50), as recommended by some authors (e.g., Quinn and Keough, 2009), identified significant correlations between SR and two descriptors: EPR and EDF.

The presence of a BEF relationship on the large spatial scale was investigated by testing the effect of different structural and functional macrofauna diversity indices on ecosystem functioning and efficiency proxies (Table 2A). The assessment of collinearity between structural and functional diversity indices with the VIF test using a threshold level of 2.5, indicated that SR exhibited a significant degree of collinearity with EPR₍₂₀₎ and EDF₍₃₀₎ (VIF = 2.7 and 2.8, respectively). This finding was expected because EPR₍₂₀₎ and EDF₍₃₀₎ are both related to the number of species (SR). Analysis of the full dataset (ATL, wM, and c-eM, see Table S5 for total macrofauna abundances) showed that macrofaunal SR was the only diversity index to be positively related to total benthic biomass, and that ES₍₅₀₎ was related to macrobenthic predator biomass (Figure 2A, Table 3). The exponential model was the one that best described these relationships (Tables 2A, 3). The relationships between other diversity indices and benthic biomass were explained by the environmental cofactors (water depth, longitude, food availability, and grain size; Table 4A).

The energy flow through the biotic ecosystem was quantified using the ratio of macrobenthic biomass to the amount of biopolymeric carbon (MBM:BPC), the ratio of macrofaunal biomass to prokaryotic biomass (MBM:TPB), the ratio of macrofaunal biomass to meiofaunal biomass (MBM: MEB), and the ratio of meiofaunal biomass to prokaryotic biomass (MEB:TPB) as proxies. On the large scale, macrobenthic biodiversity was not significantly related to any of these proxies of ecosystem efficiency (Table 2A), and most of its variability was explained by environmental covariates (Table 4A).

Macrofaunal functional diversity was expressed as trophic diversity, EDF₍₃₀₎, EPR₍₂₀₎, and BP. A BEF relationship was found only when EDF₍₃₀₎ was considered, and it was significant and exponential (Figure 2B, Tables 2A, 3). None of the other functional diversity indices used had any effect on ecosystem functioning, or else the relationships were explained by a covariate effect (Table 4A).

On the large spatial scale, BP was the only functional parameter to be positively and linearly correlated with ecosystem efficiency as measured by the MBM:BPC ratio (Tables 2A, 3). In all the other cases (Table 4A), its effect was overridden by covariate effects.

In summary, 4 of the 42 BEF relationships tested in the study indicated that macrofauna diversity significantly affected ecosystem functioning and efficiency (Table 2A). However, the collinearity between SR and EDF₍₃₀₎ suggests that their number may in fact fall to 3/42.

The relationships identified on the basin scale among macrofaunal structural and functional diversity, and ecosystem functioning and efficiency are reported in Table 2B.

Each slope system yielded different results (Table 5). The number of significant relationships was highest in the w-M basin, where macrofaunal diversity was positively related to ecosystem function and efficiency (Figure 3, Table 5). In the other slope areas, the effect of the environmental variables attenuated the BEF relationships (Table 4B). The nature of the relationships ranged from linear to exponential, according to the proxies that were applied to quantify biodiversity (Table 2B). In regard to macrofaunal functional diversity, significant and exponential relationships were found among EPR₍₂₀₎, EDF₍₃₀₎ and ecosystem functioning in the w-M and ATL area (Figures 4A–C, Table 5). In the c-eM basin, macrofauna diversity did not significantly affect ecosystem functioning and efficiency (Table 2B). Macrofaunal functional diversity only showed an exponential relationship with the MBM:MEB ratio, one of the proxies of ecosystem efficiency, and only in the wM basin (Table 5). A null relationship, that was mostly explained by the effect of environmental factors in the c-eM basin, was found in the ATL area (Table 4B). According to the present findings, the environmental variables governed all the BEF relationships in the c-eM basin and most of the relationships detected in the ATL area (Table 4B). In the wM slope area the environmental effect seemed to be slightly less marked despite the small number of significant relationships out of the number of BEF relationships assessed (Table 2B).

Continental slopes are ideal sites to investigate BEF relationships, since they account for more than 20% of total marine productivity and for a significant proportion of the organic matter exports to the seafloor (Levin and Sibuet, 2012). Slope sediments host a large quantity of marine biodiversity and are repositories of deep-sea biomass, and their benthic communities are characterized by a high spatial heterogeneity (Coll et al., 2010). The factors commonly invoked as drivers of the spatial distribution and diversity of deep-sea benthic fauna include substrate heterogeneity (Etter and Grassle, 1992), water circulation (Tyler, 2003), oxygen availability (Levin and Sibuet, 2012), productivity and microbial activity (Tyler, 2003), and food resources (Menot et al., 2010). Since the stations sampled in this study lay along a west-east axis encompassing different trophic and oceanographic conditions, environmental factors were expected to influence BEF relationships. The large spatial-scale analysis, spanning the area from the Atlantic Ocean to the Central-Eastern Mediterranean Sea, detected two significant and exponential relationships between macrofaunal diversity and ecosystem functioning. An exponential relationship between BEF has been described for various organism size classes (Mora et al., 2014); positive interspecies interactions between organisms, such as facilitation, have been suggested to sustain such relationships (Danovaro et al., 2008a). However, the diversity indices used in the present study were not all significantly related to the ecosystem functioning measures. In fact, the existence of a BEF relationship appeared to be closely linked to the diversity and ecosystem functioning measures used (Gamfeldt et al., 2014), which are often context-dependent (O'Connor and Donohue, 2013). These findings are in line with other work (Pape et al., 2013; Cusson et al., 2015; Poorter et al., 2015), where not all diversity measures correlated with ecosystem functioning.

On the large spatial scale, macrofaunal diversity—ecosystem efficiency relations were not significant. We chose our ecosystem efficiency proxies following the idea that the positive influence exerted by biodiversity on ecosystem efficiency can lead to an efficient use of available food sources and their conversion to higher biomass values (Naeem et al., 1994). Indeed, the ratio of macrobenthic biomass to biopolymeric carbon has been suggested as a proxy of ecosystem efficiency (Danovaro, 2012). This proxy has been reported to have both a positive relationship with benthic diversity (Danovaro et al., 2008a; Zeppilli et al., 2016) and no relationship at all (Leduc et al., 2013). The ratio of macrobenthic to microbial biomass and of macrofaunal to meiofauna biomass can also be employed as proxies to quantify the energy flow through an ecosystem, i.e., how efficiently it works (Cardinale et al., 2012): the higher the ratio, the more efficient the system. Similarly, the meiofaunal biomass to prokaryotic biomass ratio was used as a further proxy of ecosystem efficiency, because of the efficient carbon recycling and strong contribution to total secondary production provided by the meiofauna (De Mesel et al., 2006). Yet, all these proxies may be a gross oversimplification, because the energy flow through an ecosystem is a complex and commonly indirect flow from the smaller to the larger organisms, and is influenced by numerous biotic interactions (Piot et al., 2014), abiotic variables (Snelgrove et al., 2014; Tilman et al., 2014) and the life history of present organisms (Cusson and Bourget, 2005). This may explain the lack of significant diversity—ecosystem efficiency relationships. Moreover, we hypothesize that including a larger number of sites could lead to potentially different BEF relationships.

Because functional diversity provides a direct, mechanistic link between diversity and ecosystem functioning, numerous studies have found that functional diversity affects ecosystem functioning more strongly than taxon diversity (Reiss et al., 2009; Mineur et al., 2014). The functional diversity of organisms can be expressed in several ways, and different effects on ecosystem functioning can be identified depending on the measure employed (Naeem et al., 1995; Hiddink et al., 2009). Most of the functional diversity indices used in this study rely on the trophic diversity of macrobenthic organisms [trophic diversity index, EDF₍₃₀₎ and EPR₍₂₀₎], a key functional trait (Dunne et al., 2002; Carvalho et al., 2013). A significant and exponential BEF relationship was detected between EDF₍₃₀₎ and total benthic biomass, and it was also identified between SR and total biomass. This finding may be wholly explained by the collinearity between SR and EDF₍₃₀₎. Unexpectedly, however, while collinearity existed there was no relationship between EPR₍₂₀₎ and benthic biomass. This can probably be attributed to the fact that deposit feeders, which are the most abundant trophic group in the samples collected, may play a far from negligible role in ecosystem functioning. No relationship, but merely a weak positive trend, was detected between EPR₍₂₀₎ and predator biomass; indeed, higher numbers of predator species did not correlate with higher biomass values. Moreover, there was no correlation between predator number (ind m⁻²) and their biomass (R² = 0.03, p > 0.05); in particular, the wM slope systems were characterized by a large number of predators and a high EPR₍₂₀₎, whereas their biomass values were lower than those measured in the Atlantic slope area. Mediterranean benthic organisms have been reported to have a generally smaller size than those inhabiting the Atlantic Ocean (Sardà et al., 2004; Tecchio et al., 2011; Pape et al., 2013), and significantly lower macrobenthic biomass values have been described in the Mediterranean basin compared with the Atlantic sector (Tselepides et al., 2000; Coll et al., 2010; Baldrighi et al., 2014). The bioturbation activity of organisms can affect both the abiotic and biotic components of a system (Quéiros et al., 2015), and is one of the functional traits of benthic organisms that may sustain mutualistic interactions on the basis of BEF relationships (Loreau, 2008). On the large spatial scale, bioturbation was positively and linearly correlated with ecosystem efficiency. The biogenic activities (e.g., burrowing, reworking, and displacing of sediment particles) of benthic organisms are fundamental processes with implications for a wide range of ecosystem-related functions, such as organic matter recycling (Thurber et al., 2014). The linear relationship identified in this study indicates that all organisms contribute similarly to ecosystem efficiency (Naeem et al., 1995). This finding seems to support the notion that bioturbation facilitates organic matter recycling and uptake by higher trophic levels (Quéiros et al., 2015).

On the large spatial scale, the significant relationships we identified between macrofauna's structural and functional diversity and ecosystem functioning and efficiency were positive and linear or exponential. This result is in line with previous findings (e.g., Mora et al., 2014; Zeppilli et al., 2016), supporting the idea of a positive influence of deep-sea (macro-) faunal diversity on the functioning of ecosystem. An assumption of this work was that the steep environmental gradients characterizing the Atlantic—Central-Eastern Mediterranean transect would influence large-scale BEF relationships (Cusson et al., 2015). Since the values of BEF measures can be affected by environmental variables, it is conceivable that BEF relationships reflect the covariation of different factors, as found by the present study. Indeed, temperature, food availability, and substrate heterogeneity are among the features that determine the distribution and diversity of deep-sea benthic populations on the broad spatial scale (Levin et al., 2001). Similarly, the lack of significant associations between functional diversity (except for EDF₍₃₀₎) and ecosystem functioning found in this study may therefore be due to the influence of environmental variables on the majority of macrofaunal functional and ecosystem functioning measures (Cardinale et al., 2000; Hiddink et al., 2009). Belly and Snelgrove (2016) showed how macrofaunal functional diversity and environmental variables can equally explain the majority of variability in some key ecosystem processes. Most deep-sea communities are usually food-limited (Levin et al., 2001), and their reliance on the relatively small input of organic carbon from the surface and the particulate organic carbon flux are held to exert a strong influence on deep-sea biodiversity (Tittensor et al., 2011) and ecosystem function (Smith et al., 2008). Moreover, besides food availability, gradients in other drivers such as sediment heterogeneity, oxygen availability, and hydrodynamic regimes (Levin et al., 2001) are also capable of affecting local diversity by influencing the rates of local processes (Rex and Etter, 2010) and consequently BEF relationships.

Notably, the functional diversity indices employed did not describe ecosystem functioning more exhaustively than the traditional biodiversity indices. A possible explanation for this is that they do not encompass the full array of key macrobenthic functional traits that underpin ecosystem functioning and efficiency processes. Yet another more targeted approach, which focuses on biological traits analysis and has been proposed for the study of functional diversity (Bremner et al., 2003; Bremner, 2008; Bolam and Eggleton, 2014), also entails some limitations since our understanding of benthic functional diversity is still scant (Bolam et al., 2010). Furthermore, identification of the traits that have a direct and indirect functional role in ecosystems is highly problematic (Pakeman, 2011). The limitations are even greater where deep-sea organisms are concerned, given the poor state of current knowledge regarding their diversity and functional diversity (e.g., Smith et al., 2008; Ramirez-Llodra et al., 2010; Levin et al., 2015).

It has been suggested that BEF relationships are spatial-scale- and context-dependent (Tylianakis et al., 2008; Zeppilli et al., 2016), and that their nature is related to the system analyzed and the organisms involved (Ieno et al., 2006; Poorter et al., 2015). Analysis of the basin scale data documented a slightly different situation for each of the three slope systems. In the w-M basin, macrofaunal diversity exerted a positive effect on ecosystem function and efficiency, whereas in the other areas the effect of environmental variables explained all (c-eM basin) or nearly all (ATL) BEF relationships. The nature of the relationships identified in the Western Mediterranean basin ranged from linear to exponential, depending on the proxies that were applied to quantify biodiversity. However, irrespective of their nature, some BEF relationships were detected in the western basin but not on the large scale. According to Barnes and Hamylton (2015), the spatial scales on which biodiversity interacts with ecosystem functioning are crucial to understand the significance of BEF relationships. Also in our case, it may be speculated that some BEF relationships identified on the basin scale (the w-M basin) were not detected on the larger scale due to a masking effect exerted by environmental features. Environmental drivers, number of species, and functional diversity affect ecosystem functioning in different ways and with different strength, based on the spatial scale (Cardinale et al., 2007; Poorter et al., 2015). However, a significant relationship between macrofauna diversity and ecosystem functioning was not found in the other two areas (ATL and c-eM basin), as on the large scale. In the deep Mediterranean Sea, the west-east gradient of decreasing surface water productivity is reflected by both a decreasing amount of food reaching the seafloor (D'Ortenzio and Ribera d'Alcalà, 2009), and in a decreasing abundance of deep-sea benthic populations (Coll et al., 2010; see also Table 2). The environmental effects were strong in the c-eM basin, where most relationships were context-dependent. Indeed, environmental conditions, such as food depletion or current regime have been reported to be major factors influencing and structuring benthic populations in the Eastern Mediterranean (Tselepides et al., 2000; Króncke et al., 2003; Tecchio et al., 2011; Giovannelli et al., 2013). Furthermore, Zeppilli et al. (2016) showed that the high sea floor heterogeneity along an open slope system may shape BEF relationships, even on a very small spatial scale. In our study, we considered seven slope systems spread across a wide longitudinal range from the Atlantic Ocean to the east Mediterranean basin. The heterogeneity of the substrate was assessed through the grain size analysis, but we know that granulometric analysis cannot completely catch the seafloor morphologies even on a small spatial scale. Indeed, these seafloor structures are pivotal in shaping the BEF relations, by selecting the species responsible for the effect of ecosystem functioning (Zeppilli et al., 2016). Nonetheless, the presence and shape of BEF relationships can vary when different components (meiofauna, macrofauna, fish) are taken into account (Pusceddu et al., 2014a). This suggests that different environments (i.e., basins) may involve different functional structures of macrobenthic communities. O'Connor and Crowe (2005) concluded that different species play idiosyncratic roles, explaining why no relationship is found between species richness and ecosystem functioning in some cases. Furthermore, the inclusion of more sampling sites in the Atlantic area might lead to different BEF relationships.

The handling Editor declared a shared affiliation, though no other collaboration, with one of the authors GD, and the handling Editor states that the process met the standards of a fair and objective review. The other 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.