Unraveling the effects of environmental drivers and spatial structure on benthic species distribution patterns in Eurasian-Arctic seas (Barents, Kara and Laptev Seas)

In times of accelerating climate change, species are challenged to respond to rapidly shifting environmental settings. Yet, faunal distribution and composition are still scarcely known for remote and little explored seas, where observations are limited in number and mostly refer to local scales. Here, we present the first comprehensive study on Eurasian-Arctic macrobenthos that aims to unravel the relative influence of distinct spatial scales and environmental factors in determining their large-scale distribution and composition patterns. To consider the spatial structure of benthic distribution patterns in response to environmental forcing, we applied Moran’s eigenvector mapping (MEM) on a large dataset of 341 samples from the Barents, Kara and Laptev Seas taken between 1991 and 2014, with a total of 403 macrobenthic taxa (species or genera) that were present in ≥ 10 samples. MEM analysis revealed three spatial scales describing patterns within or beyond single seas (broad: ≥ 400 km, meso: 100–400 km, and small: ≤ 100 km). Each scale is associated with a characteristic benthic fauna and environmental drivers (broad: apparent oxygen utilization and phosphate, meso: distance-to-shoreline and temperature, small: organic carbon flux and distance-to-shoreline). Our results suggest that different environmental factors determine the variation of Eurasian-Arctic benthic community composition within the spatial scales considered and highlight the importance of considering the diverse spatial structure of species communities in marine ecosystems. This multiple-scale approach facilitates an enhanced understanding of the impact of climate-driven environmental changes that is necessary for developing appropriate management strategies for the conservation and sustainable utilization of Arctic marine systems.


Introduction
Understanding spatial patterns and temporal dynamics of biotic community structure and processes in relation to intrinsic factors and environmental forcing is a longstanding concern in ecology ( Legendre 1993). Particularly, gaining insight in the spatial distribution of species, which is not random but largely determined by intrinsic biotic properties (population demography, behavioral traits, and interspecific interactions) and external forcing (regionally heterogeneous environmental conditions) remains a challenge (Dray et al. 2012). As in all oceans, warming, acidification and shifting circulation patterns are key consequences of climate change in Arctic seas. In addition, anthropogenic activities, including fossil-fuel extraction, industrial emissions, shipping and tourism, are constantly growing and exerting mounting pressures on Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0030 0-020-02737 -9) contains supplementary material, which is available to authorized users. marine ecosystems (Doney et al. 2012;Solovyev et al. 2017). The reduction of sea ice, which is one of the most obvious effects of warming in polar regions, and the alterations of its dynamics lead to large-scale changes in the environmental conditions of Arctic seas, with cascading impacts on marine ecosystems (Macias-Fauria and Post 2018). This includes the composition and spatial distribution of benthic communities (Doney et al. 2012), the knowledge of which is of high value for environmental protection efforts (Sirenko 2001). Generally, benthic biotas are important for the overall functioning of marine ecosystems, as they decompose organic material, contribute to nutrient cycling and serve as food source for higher trophic levels (Silberberger et al. 2018). Changing environmental conditions can lead to shifts in benthic communities, such as composition and distribution, and may impact ecosystem functions (Carroll et al. 2008;CAFF 2017). Ice cover, bathymetry, water temperature, salinity, and sediment properties are considered the environmental drivers that are most influential in determining benthic biodiversity patterns and processes , while fisheries, persistent contaminants, industrial development and shipping are key anthropogenic drivers (Yasuhara et al. 2012;CAFF 2017).
Community structure varies among every existent spatial scale (broad, intermediate, and small) and is absolutely independent of organism size. Each spatial scale displays a characteristic composition of driving environmental processes and fauna. Roy et al. (2014) provide an exception to this statement in their study about the Canadian epibenthos by suggesting that there is no large regional difference in benthic community characteristics. A number of studies have revealed the importance of considering spatial scales in understanding the causes of species distributions by suggesting that ecological processes operate on multiple spatial scales and thus structure distinct aspects of species communities (Levin 1992;Bellier et al. 2007;Ingels and Vanreusel 2013;Silberberger et al. 2018;Gutt et al. 2019). However, spatial scales have rarely explicitly been considered in large-scale Arctic studies, which focused rather on providing inventories of Arctic benthic biodiversity and identifying knowledge gaps Piepenburg et al. 2011;Renaud et al. 2015). Furthermore, as sampling data for multi-scale statistical analyses are very timeconsuming and expensive, this field of research has only been scarcely investigated (Yasuhara et al. 2012). While there are several studies on local community structures in Arctic seas, few addressed explicitly multiple spatial scales in their analysis, e.g., in the Fram Strait (Budaeva et al. 2008), in the Barents Sea (Carroll et al. 2008), the Chukchi Sea (Blanchard and Feder 2014) or around the sub-Arctic Lofoten-Vesterålen islands (Silberberger et al. 2016(Silberberger et al. , 2018. A comprehensive community analysis over several Arctic seas, including the consideration of environmental drivers and a range of spatial scales, is still scarce (Roy et al. 2014).
Inventories of marine benthic fauna in the Russian Arctic have already begun at the end of the eighteenth century and have been updated regularly (e.g., Zenkevitch 1963;Sirenko and Piepenburg 1994;Sirenko 2001;Bluhm et al. 2011;Piepenburg et al. 2011;Renaud et al. 2015). Additionally, there are studies concentrating on benthic biodiversity in the different Eurasian-Arctic seas, such as the Barents (Cochrane et al. 2009;Jirkov 2013), Kara (Jørgensen et al. 1999;Deubel et al. 2003;Galkin and Vedenin 2015;Vedenin et al. 2015), Laptev (Sirenko et al. 2004;Kokarev et al. 2017) and the East Siberian Seas (Spiridonov et al. 2011). Moreover, biogeographical regionalization efforts of benthic fauna were conducted that allowed the classification of faunal species into distinct groups, such as Boreal or Atlantic subareas, and thus deliver a further evaluation of the Eurasian-Arctic marine ecosystem (Zenkevitch 1963;Spiridonov et al. 2011). However, to advance from basic regional classifications, we make use of Moran's eigenvector mapping , to reveal spatial categorizations of macrobenthic species groups and evaluate the specific impact of certain environmental drivers on faunal distribution patterns.
The specific objectives of this study were (a) detect characteristic spatial scales of benthic distribution patterns across a range of Eurasian-Arctic seas and uncover the regional distribution of the taxonomic groups in the study area, (b) identify which taxa are associated with which distinct scale, (c) assess the relative importance of environmental drivers of benthic communities at different spatial scales, and (d) assess the relative contributions of environmental factors versus spatial structure in determining the composition of benthic assemblages. This work shall add a further valuable piece of information to the global pool of existing knowledge about species-environment interdependencies in the sensitive Arctic region and help to design and implement solid science-based conservation strategies.

Study region
The area considered for our study ranges from 68° N and 16° E to 82° N and 150° E and encompasses the Barents (including the waters around Svalbard: 81.57° N 15.41° E and 68.60° N 51.00° E), Kara (81.57° N 52.30° E and 68.60° N 105.40° E) and Laptev (81.57° N 107.24° E and 68.60° N 149.46° E) Seas (Fig. 1). The majority of the Arctic continental shelves is very broad with average depths between less than 200 m in the Kara Sea (Zenkevitch 1963), 230 m in the Barents Sea (Jirkov 2013), and 533 m in the Laptev Sea (Spiridonov et al. 2011). The water column in Arctic seas usually exhibits a pronounced stratification as a result of freshwater inflow from ice melting and river discharge (CAFF 2017), particularly in the Kara and Laptev Seas, which receive 1,700 km 3 of freshwater per year from the rivers Ob, Yenisei and Lena (Vedenin et al. 2018). Between these and other Arctic rivers, like Kolyma, Pechora and Northern Dvina, large sediment loads and terrestrial organic matter are transported to the seas (Gordeev et al. 1996;Bluhm et al. 2011). Ocean currents from the Atlantic and Pacific cause mixing and exchanging processes of, e.g., nutrients, organic matter, and larvae of invertebrate species. The Atlantic feeds huge amounts of warm and salty water masses via the broad opening of the Fram Strait to the Arctic Ocean, which is supposed to affect benthic ecosystems in the study area (Rudels 2012;CAFF 2017). The Pacific Ocean supplies less warm and less salty water via the narrow Bering Strait in the Eastern Arctic, which makes it less important to the fauna in our study area (Woodgate 2018). The structure of the water column can be subdivided into three layers: an upper mixed layer which is relatively warm with temperatures of up to 12 °C and least saline, an intermediate layer with Atlantic warm and saline water and a deep Arctic layer which is colder (between 0 and − 1 °C) and more saline (between 34.9 and 35.0) (Steele et al. 2008;Rudels 2012;Woodgate 2018). The Arctic seabed consists mainly of silt and clay with the exception of drop stones providing hard substrata as habitat and ridges and plateaus which exhibit higher loads of sand ).

Biotic data
The data used in our analysis contained information on the presence of macrobenthic taxa (organisms larger than 0.5 mm) in grab samples (box corer, Petersen grab or van Veen grab, with sampling areas from 0.1 m 2 to 0.25 m 2 ) taken at a total of 341 stations in the Barents Sea + Svalbard area (145), Kara Sea (38) or Laptev Sea (158) (Fig. 1) visited in the course of 16 scientific cruises during summer (May-October) between 1991 and 2014 (Online Resource 1). The data were not taken by us, but compiled from publicly available data platforms. The prime reason why we decided to analyze presence/absence data is because they are more robust to bias introduced by varying sampling effort than quantitative data. Our biotic data rely on grab samples, which generally cannot be obtained from sandy and gravelly areas and thus are not considered in our study. Accordingly, our macrobenthic species are mainly infaunal communities as they are represented in grab samples. Sampling locations are situated approximately between 4 and 700 km off the coastline and at water depths less than 1000 m (Fig. 1). To clean and properly prepare our data for analysis, we followed the suggestions presented by Piepenburg et al. (2011). We used only those 403 macrobenthic taxa (from a total of 2,087) that (a) were identified to species or genus level and (b) occurred at least at ten sampling stations in the study area, to reduce the stochastic noise introduced by the random occurrence of 'rare' species in our data, which would bias the detection of spatial patterns. Taxonomic information was validated via the World Register of Marine Species (WoRMS; https ://www.marin espec ies.org/) to avoid the use of synonyms and outdated scientific names.

Environmental data
For our analysis, we considered 12 environmental parameters to be relevant as potential drivers of the distribution patterns of macrobenthic species in our study region  (Table 1). For each parameter, we extracted monthly average values corresponding to the place and time of sampling of the biological data. We intentionally ignored the environmental data actually measured during the respective sampling cruises, since those are mainly snapshot measurements. Moreover, for many stations there was no environmental information available from the cruises, so we used uniform datasets for all stations from publicly accessible databases. For estimating the organic carbon flux to the seafloor, using the empirical function proposed by Suess (1980), data on primary production were taken from the Ocean Productivity website of Oregon State University (https ://www.scien ce.orego nstat e.edu/ocean ) and compiled using a Carbon-based Productivity Model by Behrenfeld et al. (2005). All environmental data layers were transferred to a unified raster graphical format with the coordinate system WGS84 and rescaled using bilinear interpolation to have the same spatial dimension (67 × 552 pixels), extent (longitude between 14° E and 152° E; latitude between 66.2° N and 83° N) and latitudinal/longitudinal resolution (0.25° × 0.25°) (Online Resource 3). The reliability of data varies among variables and scales due to the distinct databases we used for data retrieval. This introduces a certain bias into the analysis, which cannot be quantified or avoided.

Multivariate spatial analysis
To determine the spatial structure in community and environmental data, a variety of multivariate modeling approaches, such as multi-scale pattern analysis (MSPA), spatial eigenfunctions or linear model of coregionalization (LMC), have been developed (Bellier et al. 2007;Jombart et al. 2009;Dray et al. 2012;Legendre and Legendre 2012). The goal of these models is to understand distributional patterns by detecting and identifying characteristic spatial scales and unraveling the relative importance of spatial distance vs. environmental forcing in explaining the spatial variation of community composition (Dray et al. 2012). Borcard and Legendre (2002) introduced the method of principal coordinates of neighbor matrices (PCNM), which in a generalized form is called Moran's eigenvector mapping (MEM; Dray et al. 2006Dray et al. , 2012. These methods are eigenfunction-based and process two types of information: binary information about whether two sampling locations are connected or not and information indicating how strong the above-mentioned connection is . We used Moran's eigenvector mapping (MEM) for the concomitant analysis of community composition, spatial structure, and environmental variables. MEM describes the different spatial scales prevailing in the distribution patterns of both the biotic communities and environmental factors. It also identifies taxa which are particularly associated with each distinct spatial scale. MEM considers implicitly the spatial distances between sampling stations in contrast to other methods which ignore this important regional feature. Moran's eigenvectors are centered, scaled, orthogonal and uncorrelated, all of which are excellent properties  (Jombart et al. 2009). The basic approach of MEM analysis is that spatial distances between sampling locations can be transformed into explanatory variables containing specific information about the relationships at different spatial scales (Peres-Neto and Legendre 2010). Creating the spatial weighting matrix (SWM), which describes the specific relationship between sampling sites, is the most crucial step of the MEM analysis. The SWM influences the accuracy of parameter estimation and the spatial patterns that can be detected with the model (Bauman et al. 2018). The SWM actually features two separate matrices: a binary connectivity matrix that indicates which sampling sites are geographically linked (1) and which are not (0), and a weighting matrix that contains information about the intensity of the connection between sites Borcard et al. 2011;Legendre and Legendre 2012).
The weighting values of the spatial links can be calculated by functions of spatial distances, least-cost links between the sampling sites or other measures (Dray et al. 2012). We tested Delaunay triangulation, Gabriel graph, relative neighborhood graph, and minimum spanning tree Borcard et al. 2011). The criteria for choosing the optimal SWM was a high level of AdjR 2 for the connectivity matrix candidates. In our case the nearest neighbor matrix best captured the connectivity between sampling sites with the highest AdjR 2 of 24.61%.
Based on the chosen SWM, we calculated the spatial eigenvalues, also called MEM variables. Since we were only interested in eigenfunctions that represented positive autocorrelations, we ran a forward selection at 999 permutations with α = 0.05 and received 49 significant MEM variables. Out of these 49 variables, 17 MEM variables displayed spatial patterns (Table 2). Correlograms were used to identify the spatial scale that each of the 17 MEM variables captures. Then, based on the similarity of their range, these variables were categorized into variable groups representing the different spatial scales that structure the community composition in the study area (Benedetti-Cecchi et al. 2010). We used forward selection to determine significant environmental variables, using a double-stopping criterion: the alpha significance level of 0.05 and the adjusted coefficient of multiple determination (AdjR 2 ) calculated using all environmental variables (Blanchet et al. 2008;Borcard et al. 2011). Afterwards, we ran a redundancy analysis (RDA) to identify and rank those variables that explained most of the variation in the community data . Forward selection and RDA were also applied to identify the associated macrobenthic fauna for each spatial scale.
Finally, we used variation partitioning to determine the shares of spatial distances and environmental factors in the explanation of spatial community patterns (Borcard et al. 1992;Peres-Neto et al. 2006;Peres-Neto and Legendre 2010). It calculates the variation in separate fractions, i.e., the variation explained by spatial distance alone, the variation explained by environmental variables alone, the variation explained by both spatial distance and environmental variables jointly, as well the fraction of unexplained variation .
Statistical analysis was conducted in R version 3.4.1 (R Development Core Team 2017), using the packages adespatial (Dray et al. 2019) and vegan (Oksanen et al. 2018).

Spatial scales and regional distribution of taxonomic groups
According to MEM analysis the distribution of macrobenthic fauna in Eurasian-Arctic seas showed distinct features at spatial scales of ≥ 400 km (broad), 100-400 km (meso), and ≤ 100 km (small) ( Table 2, Online Resource 4). MEM detected the smallest spatial distance at 54 km and the largest distance at 839 km ( Table 2).
The investigated taxa are from a broad variety of phyla, e.g., sponges, crustaceans, bivalves, polychaetes, gastropods, and echinoderms. The most frequent species, the polychaete Terebellides stroemii, was present at more than 50% of the sampling sites (57%; n = 196) (Online Resource 2). Macrobenthic fauna can be categorized into 24 taxonomic groups. In the region of Svalbard and the Barents Sea, a total of 383 taxa from 24 taxonomic groups were recorded. The Kara Sea housed 358 taxa from 23 different taxonomic groups, and the Laptev Sea a total of 291 taxa from 18 taxonomic groups ( Fig. 2; Online Resource 5). The community composition on class level was largely similar among all three regions with polychaetes, malacostracans, bivalves, bryozoans (Gymnolaemata) and gastropods being most frequent.

Macrobenthic taxa associated to spatial scales
Of the 403 macrobenthic taxa considered in our data in total, the distribution of 170 taxa showed significant structuring on at least one of the three spatial scales identified in the MEM analysis. For 50 taxa, this was evident on the ≥ 400km scale, for 90 taxa on the 100-400-km scale, and for 62 on the ≤ 100-km scale (Table 3; Online Resource 6). With a few exceptions, the composition of macrobenthic assemblages among the different scales were highly variable. The distribution of only a few taxa exhibited significant spatial patterns on more than one scale, and their ranking varied across scales. Of the 128 taxa showing significant spatial structuring at the broad or mesoscales, only 9 (7%) were associated with both scales. Of the 103 taxa that showed spatial structuring at the broad and small scales, only 6 (5.8%) were associated with both scales and of the 137 taxa showing significant spatial structuring on the meso and small scale, 12 taxa (8.8%) were associated with both scales (Online Resource 6). Only three species, the polychaetes Cirrophorus branchiatus and Lumbriclymene cylindricauda and the gymnolaemata Porella sp. displayed a significant association with all three spatial scales, suggesting that they are widely distributed in the study area.

Environmental drivers
All environmental variables used for the analysis appear as significant drivers on at least one of the detected spatial scales. On the ≥ 400-km scale, apparent oxygen utilization was identified as the major environmental factor structuring benthic assemblages (AdjR 2 of 15%), followed by near-bottom phosphate content (11.9%; Table 4). Eight of 12 environmental variables were found to be significant, together explaining a total of 33.9% (AdjR 2 cum) of community variation. On the 100-400-km scale, we found distance to shore (AdjR 2 of 3.8%) and temperature (3%) to be the environmental parameters for explaining most of the benthic community composition (Table 4). Except for oxygen and phosphate, all of the environmental parameters contributed significantly to the explanatory power of the model (AdjR 2 cum of 18.9%).
On the ≤ 100-km scale, the community variation was significantly correlated with five environmental factors, out of which organic carbon flux possessed the highest explanatory power with an AdjR 2 of 4.1%, followed by distance to shore (2.4%; Table 4). The redundancy analysis showed that the five significant environmental variables explained 9.9% (AdjR 2 cum) of the total variation in benthic community variation on this scale.

Environmental forcing vs. spatial distance
Overall, variation partitioning revealed that both spatial distances among sampling locations and environmental variables considered in our study explained part of the overall variation in the community data, ranging from 17.3% for the ≤ 100-km scale to 23.3% across all spatial scales (Fig. 3). The relative importance of spatial distance vs. environmental forcing varied among spatial scales. While across all scales, spatial distance was slightly more relevant than environment, indicated by AdjR 2 cum values of 6.7% vs. 5.3% (Fig. 3a), environmental forcing was consistently more important than spatial distance for all separate spatial scales, i.e., ≥ 400-km scale: 11.3% vs. 0.8% (Fig. 3b); 100-400-km scale: 11.7% vs. 5.2% (Fig. 3c); ≤ 100-km scale: 16% vs. 0.7% (Fig. 3d).
The variation proportions that were explained jointly by both environmental variables and spatial distance ranged from 0.6% on the small ≤ 100-km scale to 11.3% across all spatial scales (Fig. 3).

Discussion
In the Laptev Sea, despite the largest number of sampling locations (n = 158), we discern the smallest number of benthic species (n = 291) among the three investigated seas. This finding relates to the transport of Atlantic species by deep waters that do not reach the Eastern Russian Seas, due to the boundary islands of Severnaya Zemlya (Dilman 2009;Spiridonov et al. 2011). Pacific species generally play a minor role in the composition of Arctic communities, since they reach mainly the Chukchi Sea and only marginally the eastern part of the Laptev Sea (Sirenko 2001). The characteristic composition of benthic communities in the Russian Arctic is related to the bottom sediment composition, the productivity of the ecosystem, the individual location in the shelf zones and the extent of water exchange with the Atlantic or Pacific Ocean and with the adjacent seas and river runoffs (Sirenko and Piepenburg 1994;Sirenko 2009). The presence of estuarine bivalves and crustaceans prevails closer to river outfalls like the Lena Delta in the Laptev Sea or the Ob and Yenisei that flow into the Kara Sea. In greater depths All scales of 600 m up to 2500 m the communities are dominated by the occurrence of ophiurans, polychaetes, holothurians and echinoderms (Spiridonov et al. 2011). The species most associated on the ≥ 400-km scale is the bivalve Portlandia arctica, which is a very common benthic bivalve mollusk and endemic to the Arctic seas. This infaunal species prefers silty sediments and has first been described by J. E. Gray in 1824 (Gray 1824). The polychaete Leitoscoloplos mammosus is mostly structured on the intermediate 100-400-km scale and occurs in Northern Atlantic and Arctic waters. This deposit feeder has been detected by A. S. Y. Mackie in 1987 and prefers the shelf as habitat at depths of around 25 -110 m (Mackie 1987). The taxon most associated with the small scale is the polychaete Antinoella sp. This genus was first described by H. Augener in 1928 when studying the polychaetes around Spitsbergen (Augener 1928). Antinoella sp. favors muddy and sandy sediments and is widely distributed in Arctic waters such as the Canadian Labrador Sea, Icy Cape of Alaska, Spitsbergen, Bering Sea and Kara Sea (Pettibone 1993). T. stroemii, as the most abundant species in the area, has been identified as presumably not being one single species but a large complex of unrevealed cryptic species (Hutchings and Kupriyanova 2018;Lavesque et al. 2019).
The three scales that we detected in our analysis describe spatial patterns over different distances within our area of investigation. The largest distance was 839 km, suggesting that the ≥ 400-km scale embraces assemblages over more than just one sea, e.g., Barents Sea to Kara Sea or Laptev Sea. The intermediate 100-400-km scale refers to assemblage structures within a single sea. With 54 km as the smallest identified spatial distance, we suggest that the smallest ≤ 100-km scale delineates patterns covering even smaller realms within single Arctic seas. With ≥ 400 km to ≤ 100 km, the distance-range from the broadest to the smallest scale is relatively narrow in our case, while a comparable study from the sub-Arctic region could reveal regional differences on scales of several hundreds of kilometers down to 1 km (Silberberger et al. 2018). The sampling design always determines the range of spatial patterns detectable with such an analysis. We have a rather low sampling resolution, resulting in a relatively large smallest detectable scale.
The broad ≥ 400-km scale patterns that MEM analysis identified in the distribution of macrobenthic communities basically reflect biogeographical contrasts exceeding the extent of the Laptev and Kara Sea, a spatial scale which matches ecoregions or large marine ecosystems (LME) (Solovyev et al. 2017). LME can be described as regional units for conservation and management approaches of living marine resources and are characterized by unique undersea topography, current dynamics, marine productivity, and food-web interactions (Sherman 1991;Solovyev et al. 2017). LMEs are known to best address conservation issues on a regional basis. However, the finding of our broad scale may demonstrate the limitation of applying LMEs. The environmental variables that were most correlated with the ≥ 400-km scale patterns (Table 4) are proxies of benthic productivity/turnover (apparent oxygen utilization) and the strength of pelagic-benthic coupling (sediment phosphate content) (Table 4). These represent ecological differences along the longitudinal West-East gradient from sub-Arctic conditions and high-production regimes in the Barents Sea at one end to high-Arctic and low-production environments in the East Siberian Sea at the other end (Carmack and Wassmann 2006). The content of phosphate reflects sorption and desorption under different redox conditions. It is stored mostly in sediments and its release from sediment is provoked by the incapacity of the seafloor to absorb remobilized phosphate or by the high decomposition of organic matter at the sediment surface (Davenport et al. 2012;Link et al. 2013). The interconnection between bottom-near oxygen and distribution as well as activities of benthic macrofauna has already been substantiated earlier (Renaud et al. 2007;Kaskela et al. 2017).
The intermediate 100-400-km scale patterns that are superimposed on the broad-scale gradients represent spatial faunistic and environmental differences within the ecoregions along mostly latitudinal South-North gradients in distances-to-shoreline and near-bottom temperature, as indicated by the ranking of environmental variables in the explanation of community variation on this scale range (Table 4). The distance-to-shoreline variable obtains its relevance for bottom communities by the strong impact of the Siberian river runoffs (Sirenko and Piepenburg 1994;Sirenko 2009). Coastal zones are usually dynamic harsh environments with high temporal variability, which are characterized by potential seasonal hypoxia, brackish water in bays, numerous small islands off the coast and gradients in bathymetry, salinity and nutrient levels (Deubel et al. 2003). Riverine influence causes the formation of stable benthic communities that are accustomed to a freshwater environment. A vast number of our sampling points especially in the Kara and Laptev Sea are located in coastal proximity. Via freshwater discharge and the transportation of riverine species the rivers Ob and Yenisei lead to the characteristic community composition along the shelf in the Kara Sea and the Lena river in the Laptev Sea, respectively. One species that was found to be abundant in brackish waters off the Siberian coast is the crustacean Saduria entomon, which is proved by the appearance of Saduria sp. in our results for the intermediate scale (Sirenko and Piepenburg 1994) (see Online Resource 6). Bottom temperature as the second most influential abiotic driver of benthic distribution patterns on this scale is known to be one crucial factor in determining general climate and hydrographic settings. Temperature is relatively restricted in its range in the Arctic Ocean and is proved to strongly influence benthic community composition due to its effect on metabolic rates which in turn impacts growth and nutrient recycling processes (Yasuhara et al. 2012;Grebmeier et al. 2015). Grebmeier et al. (2015) suggest that macrobenthic fauna have a wide-ranging spatial distribution within the Arctic Ocean. We found that variables which are proxies for pelagic-benthic coupling are mostly driving faunal assemblages on large and small spatial scales and to a less extent on mesoscales. However, we acknowledge that this can be different, depending on the region, and in other studies pelagic-benthic coupling has been demonstrated to be most important on mesoscales (Blanchard and Feder 2014). However, the association to a spatial scale may also depend the size range of the spatial scale, which in our case is rather large.
The small-scale patterns revealed in the MEM analysis refer to structures in biotic and abiotic data on spatial scales of less than 100 km that are superimposed on the mesoscale latitudinal gradients. They basically reflect the ecological impacts of oceanographic features (such as fronts, eddies or marginal ice zones) or shelf geomorphology (such as shallows or trenches) with spatial dimensions of tens of km. In our study, such ecological impacts are primarily related to variation in general food availability for benthic communities and the environmental seafloor setting, as apparent in the relevance of environmental variables contributing to the explanation of small-scale spatial community structure (Table 4): carbon flux as a proxy of the strength of pelagic-benthic coupling and distance to shore. Organic carbon reaching the seafloor promotes high biomass, abundance and benthic biodiversity (Kędra et al. 2015). The distance-to-shoreline is due to the spatial extent of the influence of river discharge leading to an environmental setting influenced by riverine water and nutrients spread over larger areas than our smallest spatial range of around 50-100 km.
The results of the variation partitioning to determine the individual contribution of environmental variables and spatial structure clearly demonstrated the predominant role of the environment in terms of community composition (Fig. 3). At each spatial scale range identified through MEM analysis, environmental variables were significantly explaining the variation in community composition (see AdjR 2 cum values in Table 4). Similar findings were accomplished by Henry et al. (2013) and Silberberger et al. (2018). However, the considerable residuals in our models indicate that a large extent of spatial structure in the communities cannot be completely explained by the used environmental variables (Fig. 3). This indicates the influence of unknown biotic or abiotic drivers, which also affect the abundance and composition of macrobenthic communities such as food availability, bacterial abundance, predator-prey interaction, iceberg scouring, sediment stability and sedimentation rate (Włodarska-Kowalczuk and Pearson 2004;Vedenin et al. 2018). One also needs to keep in mind that, the biotic data we used have not been taken at the same time. This introduces a bias due to temporal variability and could lead to a larger amount of unexplained spatial variability. Our overall results noticeably exemplify the impact the environment and its warming-induced changing dynamic have on marine ecosystems and its inhabitants, which makes it even more substantial to investigate these physio-biological interactions further.

Conclusion and implications
Ecological knowledge on the Eurasian-Arctic seas has increased substantially over the past years; however, only few investigations have related spatial variation in benthic communities to environmental settings at various scales. So, the comprehensive dissection of spatial scales and influential environmental factors governing species composition and distribution along the Eurasian-Arctic shelf zone is new information and fosters the understanding of physio-biological dynamics. The findings of our study indicate that distribution and composition of Eurasian-Arctic benthic communities are controlled by a variety of environmental determinants acting on multiple spatial scales. Our results are specific to the Eurasian-Arctic region (Barents, Kara and Laptev Sea) and are not meant to be extrapolated to a pan-Arctic scope, due to largely varying environmental conditions across the Arctic Ocean. This kind of information is needed for current climate politics and conservation management discussions to develop suitable preservation approaches, to plan the sustainable use of Arctic resources, to apply adaptation techniques to reduce risk and vulnerability and to declare future marine protected areas (MPAs) (IPCC 2019). Marine systems, and especially remote polar seas, are generally still underrepresented in the global network of marine protected areas (MPAs) (Spiridonov et al. 2011;Solovyev et al. 2017). Ecological consequences of the accelerating climate change become continuously more severe, making it even more critical to protect certain marine habitats by identifying and assigning potential conservation priority areas (CPAs) of high ecological value (Solovyev et al. 2017).