Are Hotspots Always Hotspots? The Relationship between Diversity, Resource and Ecosystem Functions in the Arctic

The diversity-ecosystem function relationship is an important topic in ecology but has not received much attention in Arctic environments, and has rarely been tested for its stability in time. We studied the temporal variability of benthic ecosystem functioning at hotspots (sites with high benthic boundary fluxes) and coldspots (sites with lower fluxes) across two years in the Canadian Arctic. Benthic remineralisation function was measured as fluxes of oxygen, silicic acid, phosphate, nitrate and nitrite at the sediment-water interface. In addition we determined sediment pigment concentration and taxonomic and functional macrobenthic diversity. To separate temporal from spatial variability, we sampled the same nine sites from the Mackenzie Shelf to Baffin Bay during the same season (summer or fall) in 2008 and 2009. We observed that temporal variability of benthic remineralisation function at hotspots is higher than at coldspots and that taxonomic and functional macrobenthic diversity did not change significantly between years. Temporal variability of food availability (i.e., sediment surface pigment concentration) seemed higher at coldspot than at hotspot areas. Sediment chlorophyll a (Chl a) concentration, taxonomic richness, total abundance, water depth and abundance of the largest gallery-burrowing polychaete Lumbrineristetraura together explained 42% of the total variation in fluxes. Food supply proxies (i.e., sediment Chl a and depth) split hot- from coldspot stations and explained variation on the axis of temporal variability, and macrofaunal community parameters explained variation mostly along the axis separating eastern from western sites with hot- or coldspot regimes. We conclude that variability in benthic remineralisation function, food supply and diversity will react to climate change on different time scales, and that their interactive effects may hide the detection of progressive change, particularly at hotspots. Time-series of benthic functions and its related parameters should be conducted at both hot- and coldspots to produce reliable predictive models.


Introduction
A question that is central in ecology is: can diversity of communities predict ecosystem functions? This question has been feeding the experimental efforts and discussions of terrestrial and marine ecologists for some time [1,2]. However, very few studies have explicitly investigated the relationship between biodiversity and ecosystem functions at higher latitudes [3]. Ecosystem functioning, defined as the biogeochemical and biotic processes and interactions in an ecosystem, is strongly related to ecosystem services providing, e.g., wood or fish for human needs [4]. With global changes underway, such ecosystem services are threatened and, hence, efforts are increasing to define the role of biodiversity and its changes for ecosystem functions [5]. Particularly, hotspots of species richness are considered an insurance of functioning in the face of species loss [6]. Other studies have demonstrated the importance of resource availability modifying the diversity-ecosystem function relationship [7][8][9][10].
Polar ecosystems are of particular interest, because climate changes are affecting them faster and stronger than the ecosystems in other regions [11]. There is thus a need for particularly rapid assessment of how environmental changes may alter ecosystem functioning in polar latitudes. In contrast to most other oceans, more than half of the Arctic Ocean's total area consists of rather shallow continental shelves. Understanding shelf-environments is indispensable for a description of the marine Arctic ecosystems and their functioning. While recent reviews have achieved an inventory of benthic diversity of Canadian and pan-Arctic shelves [12,13], their benthic ecosystem functioning is understudied [14]. In soft-bottom environments, which dominate Arctic continental shelves, the degradation of organic matter and coupled oxygen and inorganic nutrient fluxes from the sediments back to water column is an important ecosystem function [15,16]. Benthic oxygen consumption is generally linked to the availability of food resources, which are often measured as sediment pigments [17][18][19], but other nutrient fluxes are less often integrated in ecological analyses. Knowledge on the remineralisation of other nutrients from the sediments in the Canadian Arctic has only recently been reported [20], and is not directly correlated to oxygen fluxes [21]; [22]. This underlines the need to investigate which factors influence benthic remineralisation function as a whole. In other habitats and experimental studies for example, the role of the number and identity of species for oxygen consumption and nutrient fluxes at the sediment-water interface have been clearly demonstrated [23][24][25]. Moreover, the functional diversity seems to be more important than number of taxa [26].
At sites with higher benthic remineralisation than known on average from the Canadian Arctic [20,27], which we define as hotspots, we have recently found high temporal variability of oxygen consumption. However, sites characterized by a generally low oxygen consumption (hereafter 'coldspots') showed lower temporal variability [28][29][30]. Time-series data from deep-sea sites also report important interannual variability of the downward export of organic matter [31], but changes in macrofaunal benthic community composition are observed only after several years or decades [32]. While temporal and regional variabilities in benthic oxygen uptake in the Arctic have so far been linked to environmental factors and food supply, it is not clear to which extent the benthic diversity affects benthic oxygen uptake or even multiple benthic fluxes [3,33].
For most areas of the Arctic, predictive models have to rely on few benthic data collected from different locations at different times, because Arctic time-series studies are even rarer than in other oceans [34,35]. But the use of such datasets to detect directional (progressive) change can be affected by the natural variability (stochastic change) of marine systems [34]. Climate forced environmental changes introduce additional temporal and spatial variability of ecosystem functioning. Here we assess how diversity (taxonomic and functional; composition, species number and abundance) and environmental factors affect spatial vs temporal variability of ecosystem functions at hotspots and coldspots in the Canadian Arctic. The objectives of this study were to distinguish between the temporal and spatial variation in a multivariate benthic ecosystem function in the Canadian Arctic, and to investigate the relation of the often used function proxies diversity and food supply with spatio-temporal variation. The resampling of study sites in the same season of different years for multiple remineralisation fluxes, diversity and sediment pigments was the key approach to allow for separating temporal from spatial variation in the diversity-ecosystem function analysis.
We test specifically the following hypotheses: (1) Benthic remineralisation function is significantly different between years at hotspots but not at coldspots, (2) food availability for the benthos (measured as sediment pigments) is significantly different between years at hotspots but not at coldspots, (3) Taxonomic community composition is not significantly different between years, (4) Functional community composition is not significantly different between years, and (5) Food supply explains temporal variation and macrofaunal community parameters (e.g., total abundance, richness) explain spatial variation in benthic remineralisation functioning. The results will allow evaluating whether diversity can serve as a reliable surrogate for benthic remineralisation function in the Canadian Arctic despite temporal variability in polar ecosystem processes.

Study region
The study covered the benthic ecosystems in shelf environments across the Canadian Arctic Archipelago from the eastern Mackenzie Shelf in the west to the North Water Polynya (NOW) in Northern Baffin Bay in the east ( Figure 1). These environments are characterised by strong seasonality, with the productive period being subjected to the timing of icemelt and increasing light duration as summer arrives.
The eastern Beaufort Sea and Amundsen Gulf are dominated by coastal shelves down to 600 m water depth. Pelagic primary production ranges from 30 to 70 g C m -2 yr -1 , indicating generally oligotrophic conditions [36]. Rather low primary production was also found in summer and fall 2005-2007 in the eastern Beaufort Sea, with daily production rates of 73 ± 37 mg C m −2 d −1 [37]. In the Cape Bathurst Polynya at the eastern boundary of the Amundsen Gulf, however, annual production rates are higher, reaching 90 to 175 g C m -2 yr -1 [38]. Ardyna et al. [37] reported daily primary production rates in summer of 159 ± 123 mg C m −2 d −1 , and intensive phytoplankton blooms related to ice-edge upwelling events were documented for coastal regions of the Mackenzie Shelf and Amundsen Gulf in 2008 [30,39]. Annual vertical POC fluxes of 1.6-1.8 g C m -2 yr -1 and 2.4 g C m -2 yr -1 were estimated at 200 m water depth for the Mackenzie Shelf and the Cape Bathurst Polynya, respectively [40][41][42]. Seafloor sediments are usually composed of more than 70% silt and clay [43].
The eastern North-West-Passage is marked by the opening of Lancaster Sound into Baffin Bay. From its western limitation at the Barrow Strait sill (125 m) the channel reaches a depth of more than 800 m in the sound itself. Tidal and bathymetry induced mixing of Pacific and Atlantic waters east of Barrow Strait allow for high primary production with rates of 251 ± 203 mg C m −2 d −1 [37] and an annual mean of about 60 g C m -2 yr -1 [44,45]. Vertical export can be high [46] and studies on benthic biomass and diversity report values that are among the highest known from the Arctic [47], but important gaps of data need to be filled for a more comprehensive description of the area [13, 44,48].
The North Water Polynya (NOW) is located in Baffin Bay north of Lancaster Sound. It opens each year depending on latent heat fluxes, ice-bridge formation and northerly winds over the Nares Strait between western Greenland and eastern Ellesmere Island [49,50]. Its productivity is considered to be the highest in the Arctic with primary production reaching values up to 250 g C m -2 yr -1 in the east [51] and 150 g C m -2 yr -1 across the region [36,37]. A significant amount of organic carbon is exported to water depths of > 200 m with highest values in the western polynya [52]. The seabed under the polynya reaches a depth of around 700 m and studies on benthic carbon turnover showed comparatively low rates in 1998 [53] and higher rates in 2008 [20]. Abundance of benthic fauna was highest in the NOW center, indicating the role of currents and possible advection processes in the food supply for the benthos [54].

Field sampling
Samples presented here were collected from nine sites (MD-C, AG-CW, LS-W, LS-E, NW-C, NW-E (hotspots) and MS-C, AG-CC, BB-N (coldspots), Figure 1) distributed across the study region in 2008 and 2009. To avoid confounding influence of seasons, the same sites were sampled in the same season each year. Sampling was conducted onboard the CCGS Amundsen between July and October during the Circumpolar Flaw Lead Study [55], ArcticNet expeditions in collaboration with the Canadian Healthy Ocean Network [56] and the Malina project (http://malina.obs-vlfr.fr/). Locations were chosen to study both hotspots and coldspots in the Canadian Arctic. Prior to sampling, a hotspot site was defined from previous knowledge (in order of importance, if available) as a site in an area known for high benthic oxygen fluxes, primary productivity, and/or vertical export (Mackenzie Delta plume, Cape Bathurst Polynya, Lancaster Sound, NOW). A coldspot site was defined as the opposite, outside such areas, which generally show low benthic fluxes. More information on the definition, categorisation and verification of benthic hotspots is found in Darnis et al. [20] as well as Kenchington et al. [27]. At each sampling event ('station'), an USNEL box corer was deployed for seafloor sediment collection. From each box core, three to five sub-cores of ten cm diameter and approximately 20 cm sediment depth were taken for assessing benthic remineralisation function (i.e., benthic oxygen demand and nutrient remineralisation) in shipboard microcosm incubations (Table 1). After incubation, the same sediment cores were passed through a 0.5 mm mesh sieve under slow running seawater. The sieve residues were preserved in a 4% seawater-formaldehyde solution for later analyses of species diversity and abundance under a dissection microscope. Sediment surface (first cm) of additional three sub-cores were stored in pre-weighed plastic vials and frozen immediately at -80 °C to determine food supply (i.e., sedimentary Chl a and

Benthic remineralisation function
Benthic remineralisation function was measured as the multivariate combination of benthic oxygen, nitrate, silicic acid, phosphate and nitrite fluxes at the sediment-water interface (benthic boundary fluxes). For this, shipboard incubations of sediment microcosms were run in a dark, temperaturecontrolled room (2 to 4 °C) for 24 to 48 h. Total sediment oxygen flux was determined as the decrease in oxygen concentrations in the water phase and was measured periodically (2 to 8 h intervals) with a non-invasive optical probe (Fibox 3 LCD, PreSens, Regensburg, Germany). To determine changes in nutrient concentrations, samples of the overlying water phase were taken at three times during the incubation, including the onset and end.
Oxygen and nutrient fluxes were determined as the slope of the linear regression of the oxygen and nutrient concentration on incubation time and corrected for solute concentration in the replacement water. A more detailed description of this method can be found in Link et al. [22,29].

Food supply
Food supply was measured as Chl a and phaeopigment concentrations. They were analysed fluorometrically following a modified protocol proposed by Riaux-Gobin and Klein [57] as described in Link et al. [29]. Two grams of wet substrate were incubated with 10 ml 90% Acetone (v/v) for 24 h at 4 °C, and the supernatant was measured in a Turner Design 20 fluorometer before and after acidification. Chl a and total pigment concentration (Chl a + phaeopigments) were determined. Quantities are expressed as microgram pigment per gram of dry sediment [µg g -1 ].

Macrofaunal diversity
Taxonomic diversity. Sediment residues from the sieved incubation cores were sorted under a dissection microscope in the lab to retrieve benthic organisms that were subsequently identified to the lowest possible taxonomic level and counted  (abundance, N). Taxa not identified to the species level were distinguished from other specimens (e.g. sp. 1) and classified as morpho-species. Where such consistency across the study region was not achieved (e.g., due to a lack of describable characters), specimens were grouped into the lowest common taxon (e.g., Sipuncula). Taxonomic richness is the number of taxa at each station (Tax S or S Tax ).

Functional diversity.
Consequently, species were classified into functional groups according to their traits in terms of feeding mode, body size, mobility and bioturbation influence ( Table 2, Table S1) [58][59][60]. Categories were chosen based on their presumed influence on benthic remineralisation. Species were allowed more than one trait for feeding mode. Trait information was retrieved from the best resources available [61][62][63][64][65][66][67][68]. For analyses of composition and richness, functional groups were treated in the same way as taxonomic entities. Functional group richness is the number of different categories of traits per station (S Func ).

Statistical analyses
We used a mixed-model PERMANOVA design to test for temporal and spatial differences in (a) remineralisation function (b) food supply proxies (i.e., sediment pigment concentrations), (c) taxonomic and (d) functional composition. The factors 'year' (two levels: 2008, 2009), fully crossed with 'regime' (two levels: hotspot, coldspot), 'sites' nested in 'regime' (six hotspots: MD-C, AG-CW, LS-W, LS-E, NW-C, NW-E, and three coldspots: MS-C, AG-CC, BB-N) and their interactions were tested. The resemblance matrices quantifying the between-replicate similarities in terms of all five standardized fluxes (O 2 and four nutrients) and the two sediment pigments were calculated based on Euclidean distances. Missing data points were replaced using the 'missing' function in PRIMER-E software [69]. Taxonomic and functional abundance matrices were Table 2. Functional traits.

Category Level
Feeding fourth-root transformed and their resemblance matrices were calculated based on Bray-Curtis similarity [69]. PERMANOVA pair-wise tests were run for significant sources of variation between the factors [70]. The significance level was corrected for multiple testing using the Bonferoni correction with α B = α/n, where n is the number of comparisons and α=0.05. Homogeneity of dispersion could not be tested for groups of the interaction terms 'year x site (regime)' using the PERMDISP routine due to the small sample size (n = 3) [70]. Instead, we determined average squared distances (Euclidian, for fluxes and pigments) or dissimilarities (Bray-Curtis, for taxonomic and functional composition) across sites within and between years for samples of hotspots and coldspots, respectively, using the SIMPER routine [69]. Multidimensional Scaling (MDS) plots were used to visualize the resemblance patterns.
A stepwise distance-based linear model permutation test (DistLM [71]) was performed to identify which subset of biotic and environmental variables best predicted the multivariate variation of five benthic boundary fluxes at 18 stations (all stations of 2008-2009). Ten predicting variables were entered into the model: sediment surface Chl a concentration, sediment surface phaeopigment concentration, S Tax , S Func , N, Shannon-Wiener index, abundance of the largest gallery burrower Lumbrineris tetraura (single species of its functional group) and abundance of the largest dominant tube-burrower group DFLHT, water depth and the date of ice-free conditions. Icefree conditions were determined from weekly ice charts for the western and eastern Canadian Arctic published by the Canadian Ice Service (CIS) available on http://www.ec.gc.ca/ glaces-ice/. A site was considered to be ice-free if ice concentrations below 1 prevailed for more than two consecutive weeks. To meet the linearity assumption for predictor variables, Chl a was ln transformed prior to analysis. No pair of variables was linearly correlated by r > 0.85 and hence all variables were retained for possible inclusion in the model. The stepwise routine was run employing 9999 permutations and using the AIC c selection criterion. The AIC c was devised to handle situations where the number of samples (N) is small relative to the number (v) of predictor variables (N/ v<40) [70]. Results were visualized with a distance-based redundancy analysis (dbRDA) [70].
Data will be made available in 2014 as: Link

Temporal and spatial variability of benthic remineralisation function
In general, benthic boundary fluxes were higher at hotspots than at coldspots and higher in 2008 than in 2009 ( Figure 2). This pattern was most pronounced for oxygen fluxes, whereas other nutrient fluxes showed more heterogeneous patterns. Sites of greatest benthic oxygen and nitrate uptakes, and silicic acid and phosphate releases were MD-C, NW-C and LS-W (all hotspots) (Figure 2).
The multivariate composition of benthic boundary fluxes was significantly different between hotspots and coldspots and between sites, with a significant interaction between years and the nested factor sites ( Table 3). The years 2008 and 2009 were significantly different at the sites LS-W, LS-E, NW-C, NW-E and BB-N (Table 4, Figure 3A). Variability within and between years was greater at hotspots than at coldspots ( Table 5).

Temporal and spatial variability of food supply
Sediment Chl a concentrations ranged from 0.04 to 32.44 µg g -1 with highest values at the hotspot sites MD-C and LS-W and lowest values at coldspots. Sediment phaeopigment concentrations were higher at hotspots than at coldspots (Table S2).
Significant differences in the multivariate composition of sediment pigments were found between hotspots and coldspots and between sites, with a significant interaction between years and sites ( Table 3). There were significant differences between 2008 and 2009 at sites AG-CW, MS-C and AG-CC (Table 4, Figure 3B). Variability within and between years was greater at hotspots than at coldspots (Table 5, Figure 3B).

Temporal and spatial variability of taxonomic diversity
We identified a total of 311 macrofaunal taxa in the sediments taken from the incubation cores (Table S1). Taxonomic richness (S Tax ) per core ranged from seven taxa at AG-CC (coldspot) up to 45 at LS-W (hotspot) ( Table S2). Lowest abundance was found at sites MD-C and AG-CW (both hotspots) and highest in the NOW sites (hotspot) ( Table S2).
Taxonomic composition of communities was significantly different between hotspots and coldspots and between sites, with a significant interaction between years and sites (Table 3, Figure 4A). The consecutive years 2008 and 2009 were not significantly different at any site (Table 4). Intra-annual and inter-annual dissimilarities of hotspot and coldspot communities were comparable (Table 5, Figure 4A).

Temporal and spatial variability of functional diversity
Taxa were classified into a total of 72 functional groups (Table S1). Number of functional groups (S Func ) per core ranged between six groups at AG-CC (coldspot) and 32 at LS-W (hotspot) ( Table S2).
Functional composition of communities was significantly different between hotspots and coldspots and between sites, with a significant interaction between years and the nested factor site ( Table 3). The two years 2008 and 2009 were not significantly different at any site (Table 4, Figure 4B). Dissimilarity within years was greater across coldspots than across hotspots and comparable in 2008 and 2009 (Table 5, Figure 4B).

Influence of biotic and environmental factors on the variability in fluxes
The best distance-based linear model (DistLM), explaining 42% of the overall variation in benthic boundary fluxes, was   Table 6). Sediment surface Chl a concentration contributes most to the explained variation (23.5%), followed by S Tax (9.5%), N (5.9%), water depth (4.3%) and abundance of Lumbrineris tetraura (3.9%).
Measures of functional diversity were not retained in the model. Variation of the first axis mainly separates coldspots from hotspots and pairs of the two years of each site ( Figure 5). The most important parameters contributing to the first axis of the dbRDA plot, explaining 63% of fitted flux variation, are sediment surface Chl a concentration and water depth ( Figure  5, Table 7). Benthic community parameters were most strongly correlated with the second dbRDA axis, explaining 27.9% of fitted flux variation (Table 7).

Discussion
In the Arctic, long-term prediction of the ecosystem function is often based on one-site one-year measures due to scarcity of data. Here, we focus on the different patterns of (nondirectional) temporal and spatial variability in benthic remineralisation function, functional and taxonomic diversity and resource availability, which we consider important terms of error for long-term predictions. We discuss the role of variability in hotspot vs coldspot regimes and finally, we discuss the limitations of a statistical model integrating environmental and macrofaunal parameters to explain directional temporal and spatial variation in benthic remineralisation function in the Arctic.

Hypothesis 1
Benthic remineralisation function is significantly different between years at hotspots but not at coldspots Simultaneous consideration of five benthic boundary fluxes confirmed that the magnitude of interannual variability in benthic remineralisation function differed between hotspots and coldspots. At hotspot sites the between-years differences were generally more pronounced than at coldspot sites, but the results are less consistent than we assumed. The composition of fluxes at different sites is heterogeneous due to complex interactions with environmental [22,72] and faunal [24,73] parameters. Change in benthic boundary fluxes from one year to another can be positive or negative in its direction,   Figures  2, 3A), both being hotspots. When remineralisation function as a whole is considered, it is therefore not surprising that the relative change from 2008 to 2009 across all hotspots or all coldspots is not significant, although differences between years are found for four of six hotspots. In fact, the interannual differences at the remaining two hotspots might not be detectable due to the strong within-site variation, but they show the tendency of a shift. The coldspot site BB-N was not only different in 2008 than in 2009, but also differed from other coldspot fluxes. High nitrite uptake of about 30 µmol m -2 d -1 may be the underlying mechanism. In this study, such high nitrite uptake has only been found at hotspots with other high fluxes (NW-E and LS-E) and typically indicates bacteriamediated anaerobic degradation of nitrogen derivates such as ammonium oxidation (anammox). Anammox is considered a common process in deep Arctic cold-water environments, and the nitrite uptake rates measured in 2008 might be a lag response to the degradation of intensive organic matter pulses fuelling abundance and degradation by anammox bacteria [74]. Site-dependent changes in nutrient fluxes were also found in a shallow-water environment by Thouzeau et al. [75], who interpreted interannual differences of biogeochemical fluxes as being mostly affected by differences in environmental parameters and organic matter deposition, including indirect effects of macrofauna and macroalgae, depending on the sites. Similar factors could have affected the changes at our hotspot sites.
One of the rare time-series studies including the measurement of total sediment oxygen fluxes in deeper water was conducted in the abyssal north Pacific for eight years [76]. Although the export of organic matter varied between years, oxygen consumption remained fairly stable. Smith et al. [76] interpreted this discrepancy in the rather oligotrophic environment they investigated as a capacity of the benthic fauna to endure food deficiency over a limited time perioduntil a new food pulse arrives. Such an explanation would also fit to oxygen flux patterns of coldspots in our study. In a less oligotrophic environment, Lepore et al. [77] compared sediment oxygen fluxes (expressed as benthic carbon respiration) between years in the Chukchi Sea in 2002 and 2004. Summed over several sites, they found that benthic carbon respiration did not change much, neither on the shelf nor on the slope over the years, despite great changes in vertical carbon export. One possible explanation for this finding is a time lag in the response of the sediment communities. However, across-site patchiness in benthic carbon respiration was high and may have masked temporal variations at particular sites -as it was the case across hotspots and coldspots in our study. This stresses the importance to separate spatial from temporal variability if we want to understand dynamics of ecosystem functions. A few sites in our study have been sampled prior to 2008. In 2004, Renaud et al. [78] reported less than half the oxygen demand (5.65 mmol m -2 d -1 ) of our values at hotspot site MD-C but similar values at hotspot site AG-CW (2.12 mmol m -2 d -1 ). At hotspot site NW-C, the oxygen fluxes we measured in 2008 were twice as high as those measured in 1998 (4.3 mmol m -2 d -1 [53]), but they were largely the same in 2009. Longer time-series measurements are necessary to draw a firm conclusion, but these results may indicate that a progressive change is already happening in the Mackenzie Delta region, whereas there is stochastic variability in the NOW region. We found that differences between sites can be in the same order of magnitude as variations across years ( Figure 3A). However, the significant difference in remineralisation function between hotspots and coldspots statistically confirms our a priori categorisation of sites. Although we cannot separate a location from dispersion effect [70] in the data, our results clearly showed that benthic remineralisation function is more variable between years at hotspot sites than at coldspot sites. This means that quantifying progressive directional changes in ecosystem functions at hotspot sites requires a long-term data series [34], whereas coldspot sites might require less regular sampling to detect changes.

Hypothesis 2
Food availability is significantly different between years at hotspots but not at coldspots Only one of six hotspot sites but two of three coldspot sites changed significantly over the years, but our results showed that the magnitude of interannual variations in food supply to the benthos, as approximated by sediment pigment concentrations, differed between hotspots and coldspots and depends on the considered site. Also, hotspots were different from coldspots. How can we explain the unexpectedly rare interannual changes at hotspots and more common changes at coldspots? Here, we analysed concentration of rather labile Chl a and stable, more degraded phaeopigments in the sediments simultaneously. Phaeopigments can accumulate with degradation of matter and are therefore not necessarily an indicator of recent organic matter input, which we want to detect when looking at annual variability. But since Chl a is often rapidly degraded to phaeopigments [79], the latter allow the detection of food input even if it had been of lower quality or earlier in the year. Another indicator of food at the seafloor is the ratio of Chl a over phaeopigment concentration (e.g. [80]). But as a quality ratio, this measure does not allow the comparison of the quantity of pigments at different sites. Since the quality of food is often related to the activity of benthic communities [29,79], we could consider only using sediment Chl a concentration, which should be representative of fresh matter input having occurred in the same season. However, statistical analyses using only sediment Chl a did not yield very different results, and only hotspot MS-C differed significantly between years (data not shown). Furthermore, one important issue that could make changes statistically less detectable is the small-scale spatial (within-site) variation in our sediment pigment data [70]. Furthermore, dispersion was greater across hotspot sites than across coldspot sites, impeding the detection of changes in the plot. In fact, the results from the SIMPER analysis clearly showed that hotspots were more different from one year to another than coldspots, but also that the difference between years for coldspots was large compared to the variability between coldspots of the same year.
Interannual changes of sediment pigments were found at several depths in the Fram Strait at the HAUSGARTEN site in Norway [35], which could be related to a decrease in the vertical flux of phytodetrital matter. In our study area, interannual variability of vertical flux patterns is only known for part of the study period and study area. In the southeastern Beaufort Sea, Sallon et al. [81] found two-day vertical fluxes at  Proportion of variance in benthic boundary fluxes explained by the 4 axis retained in the DistLM (Table 6) hotspots MS-C and AG-CW in 2008 that were similar to those reported in 2004 [82]. Both were higher than those determined in the 1980s [42]. Recent results from 2009 showed, however, that vertical fluxes at theses sites were lower than in 2008 [22,83]. Due to the generally tight pelagic-benthic coupling in the Arctic, such interannual variability in vertical fluxes in the southeastern Beaufort Sea (see also 84 for a more regional approach) could lead to interannual sediment pigment changes, even at coldspots. In the NOW, at the hotspots NW-C and NW-E, Hargrave et al. [52] showed an about two-fold increase of vertical flux from 1998 to 1999. In the LS-W region, strong interannual variability of processes determining pelagicbenthic coupling, including the release of ice algae, have been observed between 1984 and 1992 [44]. However, in that period the change seemed stochastic, and contrarily we may have sampled in two years of similar vertical export in 2008-2009. Although the different studies on all the subregions show interannual variability, we do not know how or whether vertical fluxes changed between 2008 and 2009, and it is therefore difficult to explain the lack of signal for sediment pigments using vertical flux patterns. In the area of LS-E or BB-N, no data on vertical fluxes has been published to our knowledge, thus that we can only infer variability from primary production and phytoplankton Chl a biomass patterns. Significant interannual differences between 2005 and 2007 (confounded with season late summer, early fall and fall) in primary production have not been found for the Beaufort Sea, Archipelago (including Lancaster Sound) or NOW [37]. However, phytoplankton Chl a biomass at individually resampled sites AG-CW (408), AG-CC (405), NW-C (108) and NW-E (115) seemed different in 2005 than in 2006 and 2007. But at water depths of more than 200 m, several processes influence the actual export of primary production to the seafloor. Temporal and spatial variability in these pelagic processes may nullify or even invert an increase in primary production to a decrease in vertical export, as has been reported at some sites in the southeastern Beaufort Sea in 2008 [28,81]. Such discrepancy is particularly evident for site  AG-CC, which Ardyna et al. [37] identified as a hotspot of primary production, but which is consistently a coldspot in benthic parameters. Lateral advection has been proposed as mechanism underlying discrepancies between surface primary production and benthic activity [77]. However, currents in the central Amundsen Gulf region are generally weak, thus that pelagic degradation is probably more important [28,84]. Extrapolating temporal and spatial changes in benthic food supply from primary production patterns should be treated with extreme caution. Benthic degradation of organic matter could theoretically also be the reason why coldspots with consistently lower sediment pigment concentration demonstrate more changes than hotspots: Assuming that hotspots generally host higher abundance of organisms due to generally higher food supply (e.g. [85]), such higher density could consume a peak of arriving food more rapidly than a site with lower abundance. In our study, surface deposit feeders were more abundant at hotspot sites (average of 120 per site) than at coldspot sites or sites with significant differences in food supply between years (average of 10 or 12 per site, respectively) (data not shown). Rapid consumption of organic material has been demonstrated for seafloor communities of different depths [86] and Arctic locations [29,79]. If sediment food supply is measured months after the input peak, the measured sediment pigments could represent the 'leftovers of the feast', which should vary less when abundance is high. In this case, food supply may work better as proxy for benthic functions, where rapid deposit feeding groups are absent. It is also interesting to note, that average similarity of functional community composition between years is lower at sites with significant differences in food supply between 2008 and 2009 (Table 4). Overall, the results indicate that functional community composition may bias food supply proxies. Hence, when explaining changes in benthic function, the interaction of community composition on measured food supply (i.e. the remaining stock) could influence the results.

Hypothesis 3
Taxonomic community composition is not significantly different between years Macrobenthic communities at hotspots or coldspots in the Canadian Arctic did not change significantly in taxonomic composition from 2008 to 2009. Instead, the community composition is generally different at hotspots than that at coldspots, and sites are distinct from each other. The total number of 331 macrobenthic taxa found in our samples may seem low compared to the overall number of 992 taxa reported from the entire Canadian Arctic [87]. Piepenburg et al. [13] performed a rarefaction analysis based on molluscs, arthropods, echinoderms and annelids for different regions in the Arctic. They showed that 19 sampling events would yield an average of 274 observed taxa in the Amundsen Gulf region and 205 on the Beaufort Shelf. Moreover, based on fewer stations, totals of 86 and 204 taxa have been reported for Lancaster Sound and northern Baffin Bay, respectively [13]. Considering species overlap between the different regions and a total of 18 sampling events, we can assume to have taken and analysed a typical number of species for the sampling size. The number of taxa we found at sites in the NOW region was also lower than reported by Lalande [54]. But due to the different sampling approaches applied (here: 3 incubation cores of 80 cm 2 each vs large corers covering areas of usually 1250 cm 2 ), our results are not directly comparable, and the numbers most likely only represent differences in total sampled surface. Two data distribution patterns in the MDS plot ( Figure  4A) are interesting to note: First, with the exception of site AG-CW 2008, hotspot sites are well separated from coldspot sites. This observation supports the notion that species composition is largely dependent on regimes of food supply [88], since primary production regimes were part of our hotspot definition. Second, sites are clustered according to their depth within the hotspot or coldspot groups (which is much less the case for benthic fluxes) ( Figure 4A). This indicates depth zonation of benthic infauna communities on a shelf gradient, which has previously been reported [43,89], but could not be identified by Cusson et al. [12] in other regions of the Canadian Arctic. These results of within-regime gradients stress the importance that environmental factors other than depth create communities of different diversity and composition.
In contrast to spatio-temporal patterns in benthic remineralisation function and food supply, the significant interaction between the factors year and site in our data results only from statistically different pairs of different sites across years. In support of this, similar variability of communities across sites in 2008 and 2009 is confirmed through the SIMPER analysis. Polar environments host species with longer life spans, and community turnover rates should therefore be slower. In an Arctic Fjord, Kedra et al. [63] reported a decadal shift in communities in two of the fjord regions influenced by an intrusion of warm water, but no changes in a third region with more stable hydrographic conditions. This supports the notion that community changes are rather long-term reactions to environmental forcing. Similarly, long-term time series of macrobenthic infauna in the northern Bering Sea show progressive change in community composition since 2000 on a decadal scale, although the change is not yet statistically detectable [32]. At sites in Frobisher Bay (Canadian Arctic), however, Cusson et al. [12] detected monthly variability in community composition. Considering that those communities became similar over time again, the effect of lacking replicate data information could have reduced both within-month and within-site variability and thus produced an artificial signal. The overall low within-site variability of Frobisher Bay sites compared to other sites in the Canadian Arctic supports this explanation. Interannual variations in meiofaunal abundance were found during time-series measurements at the HAUSGARTEN from 2000-2004 [89]. However, these changes were less clear when community composition was considered. Overall, the results support our hypothesis and spatial variability in community composition seems more pronounced than interannual variations. However, without separating the temporal from the spatial component of variability, we will miss detecting long-term shifts in community composition, such as the shift reported by Grebmeier [32] in the northern Bering Sea.

Hypothesis 4 Functional community composition is not significantly different between years
Our results support the hypothesis that macrobenthic communities at hotspots or coldspots in the Canadian Arctic did not change significantly in functional composition from 2008 to 2009. Instead, the functional community composition is generally different at hotspots than that at coldspots, and sites are distinct from each other. The number of 72 functional groups reported in our study may seem very high, but if functional groups are to be related to ecosystem functions, it is crucial to include as many different trait categories as known to influence those functions [60]. We therefore think that such fine-scale categorisation is adequate. Functional richness (number of functional groups) was largely similar to the number of taxa at sites with only few taxa, but was clearly lower at sites with many taxa (Table S2). This is a typical phenomenon, since in small-scale studies, the probability to find species redundant in their functions increases with the total number of species encountered [90]. Functional diversity can follow changes in taxonomic diversity -or not [91,92]. If functional redundancy is present (as for samples where taxonomic diversity is higher than functional diversity), a reduction in the number of species will not necessarily decrease the number of functional groups. This implies less variability, and the absence of changes in functional composition over the years matches well with temporal patterns in taxonomic composition.
When comparing results of variation in functional composition and taxonomic composition, we again find a welldefined separation of hotspot and coldspot sites. However, sites are generally more similar to each other when functional diversity is analysed. Using functional diversity in comparative studies has received increasing attention over the last decades [2]. An important advantage of this approach is its lower susceptibility to misclassifications, and Cochrane et al. [93] have recently demonstrated its utility to describe ecologically distinct regions in the Barents Sea. Although the general power to identify spatial and temporal differences in community composition is reduced when the functional (opposed to taxonomic) approach is used, functional diversity is easier to use and sufficient to describe the ecological role of different communities.

Hypothesis 5
Food supply explains temporal variation and macrofaunal community parameters explain spatial variation in benthic remineralisation function The results of this study confirm our hypothesis that the temporal variability of benthic remineralisation function is most affected by sediment pigments while its spatial variation is largely determined by diversity patterns. However, water depth was an additional important factor explaining spatial patterns, both across sites and between hotspots and coldspots. We demonstrated how depth played a role in explaining the resemblance patterns in taxonomic community composition (see section Hypothesis 3). Moreover, water depth has been used as a rough estimate for food supply to the benthos [94], as also indicated by its similar relationship to the variation axes as that of sediment Chl a concentrations in our data. The importance of an easy-to-determine and steady-state variable for ecosystem functions can have stabilizing effects in predicting models, but also calls for caution considering its low explicative power and when ecosystems of different depths and regimes are compared [34]. In our data, sediment Chl a concentrations are strongly related to the variation axis across which hot-and coldspots are separated, and along which temporal variability of sites is spread, despite the lack of temporal variability in the variable itself. The retention of Chl a but not phaeopigments in our model underpins the notion that fresh food supply rather than general food supply is most important in determining benthic remineralisation function [22,29,79].
About 20% of the variation in benthic function is explained by different measures of diversity (richness of taxa, individual abundance, and the abundance of L. tetraura). These three variables mostly explain differences between hotspots and coldspots and the eastern and western Canadian Arctic, but also temporal variation at hotspots AG-CW and NW-E. The better explicative power of taxonomic vs functional group richness may be due to several mechanisms. When relating functional diversity to ecosystem functions, species traits relevant for the considered ecosystem functions need to be included [60]. Bioirrigation has been shown to influence biogeochemical fluxes across the sediment-water interface [73,95,96]. The lack of information for species in our study did not allow including this trait as a category, but including traits more directly linked to functions, in our case e.g. including bioirrigation may improve the explanatory power of functional type models for ecosystem functions. Further, resource availability, measured as Chl a in our study, may affect the diversity-ecosystem function relationship. At large spatial scales across the marine Arctic, the number of benthic species generally increases with primary production, if the effect of salinity is removed from the model [97]. If functional group richness is more related to sediment Chl a concentration, then taxonomic richness is chosen as a more complementary variable to the already used Chl a in the model. Moreover, interactions between functional and species richness in their effect on functions have been found, making species richness more important, if functional richness is low [10]. Another possible cause for taxonomic richness being more important is the use of multiple processes defining our function: Using multiple processes decreases the chance that several species are redundant in effecting the functions [98].
Our statistical approach of a predictive linear model does not allow testing the combined effect of food supply and whole community composition on multiple benthic fluxes. But total abundance and the abundance of a gallery-burrowing polychaete species as significant predictors for multiple benthic functions stress the importance of community composition additionally to mere species numbers. In fact, the density of fauna [99], identity or functional traits of species [24,25] and the number of burrows at the seafloor [73] have been related to benthic boundary fluxes before.
More than half of the variation in benthic boundary fluxes in our data could not be explained by our best model. As already mentioned above, more exact measures of community composition or trait composition could explain more variation. Two other important benthic community components are also lacking in our study: meiobenthic and bacterial abundance. The role of meiobenthos for biogeochemical cycles is less studied, but there is increasing evidence that its density is related to benthic remineralisation function [25,100,101]. The role of bacteria is known to be relevant for organic matter degradation processes [102]. Different types of bacteria can influence the particular process of degradation (e.g. anammox), can react rapidly to matter input and interact with macrofaunal matter degradation processes [74,103,104]. Although evidence is increasing that macrofauna drives the variability in bacteriamediated degradation [104,105], future studies should integrate meiofaunal and microbenthic organisms to gain a better understanding of mechanisms regulating temporal and spatial variability of benthic remineralisation function.
Benthic ecosystem function in the Canadian Arctic increases mostly with fresh food supply, species-rich and functional diverse communities. With climate change, the quantity andeven more importantly -the quality of food supply to polar benthos are changing [14,28,30]. While we currently still know only little about how benthic communities will react to these changes [32], we could show here that diversity changes will also have an impact on the quality of benthic ecosystem functions. As recently indicated in a metaanalysis on studies from the aquatic environment [8], our study confirms the comparable importance of abiotic factors and diversity for ecosystem functions in a natural marine environment.

Conclusion
Great efforts are underway to estimate the impact of climate change on polar ecosystems. Due to the very limited number of real benthic time-series measurements in the Arctic, ecosystem models often rely on data obtained from different sites at different times of the year. One important question under these circumstances is: Are known benthic diversity hotspots in the Canadian Arctic also hotspots in ecosystem function? We conclude from our findings: Yes, with regard to a general spatial comparison, but No, with regard to long-term predictions, which may increase variability of resource availability due to climate-related ecological change. In our study, we have demonstrated that the influence of diversity on multiple benthic ecosystem functions is complementary or dependent on the availability of resources. The mechanisms controlling temporal variability of factors explaining benthic ecosystem function vary even on a within-region spatial scale. We have also shown that the similarity in taxonomic and functional diversity patterns indicate little insurance against climate induced species loss in the Canadian Arctic. Even models that include several variables (steady-state ones but also those varying on short and long time scales) explained only half of the variability in multiple benthic ecosystem function in the Canadian Arctic. Defining the role of the functional identity of particular organisms in benthic biogeochemical cycles should help to better predict benthic remineralisation function. Our findings also strongly suggest that for reliable predictions of how ecosystem functioning in Arctic shelf habitats will change in the future and how close we are to tipping points, it is necessary to establish time-series sites at hotspots and coldspots where multiple function measures are monitored, in order to distinguish natural oscillations from directional change.