Patterns in Benthic Microbial Community Structure Across Environmental Gradients in the Beaufort Sea Shelf and Slope

The paradigm of tight pelagic-benthic coupling in the Arctic suggests that current and future fluctuations in sea ice, primary production, and riverine input resulting from global climate change will have major impacts on benthic ecosystems. To understand how these changes will affect benthic ecosystem function, we must characterize diversity, spatial distribution, and community composition for all faunal components. Bacteria and archaea link the biotic and abiotic realms, playing important roles in organic matter (OM) decomposition, biogeochemical cycling, and contaminant degradation, yet sediment microbial communities have rarely been examined in the North American Arctic. Shifts in microbial community structure and composition occur with shifts in OM inputs and contaminant exposure, with implications for shifts in ecological function. Furthermore, the characterization of benthic microbial communities provides a foundation from which to build focused experimental research. We assessed diversity and community structure of benthic prokaryotes in the upper 1 cm of sediments in the southern Beaufort Sea (United States and Canada), and investigated environmental correlates of prokaryotic community structure over a broad spatial scale (spanning 1,229 km) at depths ranging from 17 to 1,200 m. Based on hierarchical clustering, we identified four prokaryotic assemblages from the 85 samples analyzed. Two were largely delineated by the markedly different environmental conditions in shallow shelf vs. upper continental slope sediments. A third assemblage was mainly comprised of operational taxonomic units (OTUs) shared between the shallow shelf and upper slope assemblages. The fourth assemblage corresponded to sediments receiving heavier OM loading, likely resulting in a shallower anoxic layer. These sites may also harbor microbial mats and/or methane seeps. Substructure within these assemblages generally reflected turnover along a longitudinal gradient, which may be related to the quantity and composition of OM deposited to the seafloor; bathymetry and the Mackenzie River were the two major factors influencing prokaryote distribution on this scale. In a broader geographical context, differences in prokaryotic community structure between the Beaufort Sea and Norwegian Arctic suggest that benthic microbes may reflect regional differences in the hydrography, biogeochemistry, and bathymetry of Arctic shelf systems.


INTRODUCTION
The Arctic marine ecosystem is undergoing pronounced changes due to climbing atmospheric temperatures, occurring two to three times faster than the global average (ACIA, 2005;Walsh et al., 2011). Resulting shifts in sea-ice cover, primary production, and riverine input will likely affect the quality and quantity of organic material (OM) deposited to the seafloor and the tight pelagic-benthic coupling characteristic of Arctic ecosystems may emphasize the effects of these changes (Grebmeier, 2012;Kortsch et al., 2015). Impacts of climate change on Arctic benthic eukaryotes have varied among size classes, but benthic prokaryotic communities have not been described well enough to monitor change (Nelson et al., 2014;Kêdra et al., 2015;Renaud et al., 2019).
Benthic bacteria and archaea perform diverse metabolic functions that mediate biogeochemical processes, including degradation and early diagenesis of organic matter in marine surface sediments (Deming and Baross, 1993). Prokaryotes provide a link between the abiotic and biotic realms, such that prokaryotic community structure reflects environmental gradients in, for example, OM deposition. Studies conducted in the Norwegian Arctic suggest that prokaryotic community structure shifts with varying quality and quantity of OM inputs, with possible consequences for ecosystem function (Hoffmann et al., 2017;Braeckman et al., 2018). Certain taxonomic groups have been strongly correlated with environmental parameters, such as Chl-a and bathymetry, measured along natural environmental gradients (Bienhold et al., 2012;Jacob et al., 2013). Benthic prokaryotes are also key players in both the production and degradation of greenhouse gases such as methane and nitrous oxide and represent a first line of defense in remediating contamination by oil and other petroleum-based products (Boetius et al., 2000;Head et al., 2006;Casciotti and Buchwald, 2012). In some cases, these important functions are attributable to specific taxonomic groups (Kostka et al., 2011;Ruff et al., 2015;Otte et al., 2019). Thus, characterization of benthic prokaryotic community structure may yield valuable insights into ecosystem function in Arctic marine sediments, particularly those permeated with methane and subjected to active mineral resource exploration such as the Beaufort Sea sites examined here BOEM, 2019).
The Beaufort Sea is an Arctic marginal sea that crosses a border between the United States and Canada, stretching from Point Barrow (Alaska, United States) to Banks Island (Northwest Territories, Canada). Prior studies of benthic prokaryotes in the southern Beaufort Sea have focused on subsurface sediments (20-70 cm) to assess the influence of either methane or mud volcanoes on bacterial community structure (Hamdan et al., 2013;Kirchman et al., 2014;Treude et al., 2014;Lee et al., 2018). Communities in these deeper sediments typically differ substantially from the active surface layer (Inagaki et al., 2003;Fry et al., 2006;Roussel et al., 2009;Durbin and Teske, 2011;Jorgensen et al., 2012;Bienhold et al., 2016). Thus, there is no established baseline for benthic prokaryotes in southern Beaufort Sea surface sediments to date, and only limited data from deeper sediment horizons. Most information on Arctic benthic prokaryotes comes from the Norwegian Arctic, largely from deep-sea or fjord habitats, which differ from the Beaufort Sea in terms of primary production, terrestrial/riverine input, and bathymetry (Ravenschlag et al., 2001;Holmes et al., 2002;Carmack and Wassmann, 2006;Arnosti, 2008;Li et al., 2009;Jacob et al., 2013;Buttigieg and Ramette, 2015;Bienhold et al., 2016;Zeng et al., 2017).
Across the southern Beaufort Sea, summer depth-stratified water-mass structure includes a fresher surface Polar Mixed Layer (PML; ∼0-50 m), the Pacific-influenced Arctic Halocline Layer (AHL; ∼50 -200 m), and warmer, more saline Atlantic Water (AW) at depths greater than 200 m (Lansard et al., 2012;Miquel et al., 2015;Smoot and Hopcroft, 2017). Eastern Beaufort Sea shelf (< 100 m water depth) and slope habitats are particularly influenced by riverine input, including the vast Mackenzie River which is the fourth largest in the Arctic in terms of freshwater discharge (330 km 3 /year), and the largest in terms of sediment transport (124 × 10 6 t/year) (Macdonald et al., 1998;Holmes et al., 2002;Rachold et al., 2004). The Mackenzie River primarily influences the water and sediment characteristics of the Canadian sector of the shelf and slope, and to a lesser extent the Amundsen Gulf and easternmost section of the Alaskan sector (Figure 1). The Mackenzie shelf functions as a vast estuary influenced by freshwater runoff, receiving inputs of both terrestrial and marine sources of organic matter (Carmack and Macdonald, 2002;Magen et al., 2010). The Amundsen Gulf, which connects the Canadian Arctic Archipelago with the southeast Beaufort Sea, is characterized by several peripheral bays, straights, and inlets and ringed by a narrow shelf surrounding a central basin (Stokes et al., 2006;Forest et al., 2010;Gamboa et al., 2017). Water mass properties in the Alaskan Beaufort Sea are influenced by warmer, more nutrient-rich waters entering from the Chukchi Sea and forming the Beaufort Shelf Jet, which flows due east from Barrow Canyon (Pickart, 2004;Weingartner et al., 2005;Dunton et al., 2006). These longitudinal environmental gradients from Point Barrow to Banks Island have been reflected in the benthic distributions of OM and certain epi-/infaunal groups (Naidu et al., 2000;Magen et al., 2010;Goñi et al., 2013;Bell et al., 2016). This benthic habitat may also be heavily influenced by methane derived from organic matter burial, mud volcanoes, subsurface permafrost, and/or gas hydrates (Paull et al., 2015;Brothers et al., 2016;Coffin et al., 2017;Lee et al., 2018;Retelletti Brogi et al., 2019).
We sought to describe prokaryote communities across this broad, heterogeneous benthic habitat, specifically by (1) assessing the diversity, structure, and composition of archaeal and bacterial communities in the upper 1 cm of sediments, and (2) investigating environmental correlates of prokaryote community structure over a broad area of the southern Beaufort Sea shelf and slope. Expanding knowledge of functional taxonomic groups of marine prokaryotes provides insight into this unstudied microbial habitat, and a framework for developing more targeted experimental studies in this dynamic polar environment (Fuhrman, 2009;Bier et al., 2015;Fuhrman et al., 2015). FIGURE 1 | Prokaryotic cluster assemblages and spatial distribution in Beaufort Sea sediments. This figure shows the dendrogram of the four assemblages revealed via hierarchical clustering above the study area map which exhibits sample locations colored by assemblage. The map was created using ArcGIS online with the following layers: Northwest Territories, Esri, Garmin, FAO, NOAA, USGS, EPA, NRCan, and Parks Canada.

Sample Collection
Our study area covers a broad section of the continental shelf and upper slope of the Alaskan and Canadian Beaufort Sea, from areas offshore of the Colville River in the west to Banks Island in the east, and into the Amundsen Gulf (Figure 1). A total of 85 sediment samples were collected opportunistically at water depths between 17 and 1,200 m from 2012 to 2014 as part of two international collaborative field programs: the US-Transboundary Fish and Lower Trophic Communities Project (USTB) and the Canadian Beaufort Regional Environmental Assessment Project (BREA). Fifteen sediment samples, 8 in 2012 and 7 in 2013, were collected using a 5 cm diameter sub-core from the upper 1 cm layer of a 0.25 m 2 box core. In 2014, 70 sediment samples were collected using a 60-cc, 2.5 cm diameter, sterilized syringe from the upper 1 cm surface layer of a double van Veen grab (0.1 m 2 , USTB) or a 0.25 m 2 box core (BREA). As preliminary analyses indicated that neither year or sampling gear were significant factors affecting prokaryotic community structure here, all samples were used in this study. The samples were immediately frozen at −20 • C following collection aboard respective vessels, transported to the University of Alaska Fairbanks (UAF) Institute of Arctic Biology Genomics Core Laboratory, and then stored at −80 • C. Environmental variables such as depth, bottom water temperature, and salinity were recorded from CTD profiles conducted at corresponding sampling locations Niemi et al., 2015;Smoot and Hopcroft, 2017).

Sediment Characteristics
We quantified several common indicators of quality, quantity, and origin of sediment OM including chlorophyll-a concentration (Chl-a), phaeopigment concentration (Phaeo), total organic carbon (TOC), total nitrogen (TN), ratio of carbon to nitrogen (C:N), and bulk sediment δ 13 C and δ 15 N. All parameters were measured using sediment samples taken from the upper 1 cm of the same box core or grab as the genomic samples from each sampling station. Chl-a and phaeopigment concentrations were measured fluorometrically (Arar and Collins, 1997). Briefly, samples were suspended in 5 ml 100% acetone, sonicated in an ice-water bath for 10 min, and allowed to extract overnight at −20 • C. Samples were then centrifuged to remove sediment, and transferred to a clean test tube. Chl-a concentration was determined using a TD-700 fluorometer (Turner Designs). After recording fluorescence values, samples were acidified with HCl, and fluorescence readings were taken of the acidified samples to produce phaeopigment values. A standard curve produced using commercially available Chl-a standard was used to convert fluorescence readings into concentrations.
Stable isotope and elemental analysis were performed on freeze-dried sediment samples at the Alaska Stable Isotope Facility (ASIF). Data for USTB samples were generated and previously reported by Bell et al. (2016); additional samples were analyzed here using the same protocol. Stable isotope data, δ 13 C and δ 15 N, were generated using a ThermoFinnigan DeltaVPlus isotope ratio mass spectrometer with Pee Dee Belemite (PDB) and atmospheric nitrogen as standards for carbon and nitrogen, respectively. Percent carbon and nitrogen content was obtained using a Costech ESC 4010 elemental analyzer, and used to calculate TOC, TN, and C:N ratios. Sediment porosity ( ), i.e., volume of water within sediments (V w ), was calculated as an indicator of sediment permeability/O 2 penetration using Eq. 1 (Bennett and Lambert, 1971;Ullman and Aller, 1982). Sediment wet (W sw ) and dry (W sd ) weights were recorded before and after freeze-drying, and used to calculate the mass of water (Ww). Sediment density (ρ s ) was held constant at 2.50 g/cm 3 , based on Coffin et al. (2017). Seawater density (ρ w ) was calculated based on bottom-water temperatures and salinities recorded at sampled sites.

Microbial Community Analyses
16S ribosomal (rRNA) gene amplicon sequencing was conducted on freeze-dried sediments in order to assess the diversity and taxonomic composition of prokaryotes. Total genomic DNA was extracted from sediment samples using the Qiagen PowerSoil kit, and revised forward (515FB) and reverse (806RB) primers from the Earth Microbiome Project (EMP) were used to amplify the V4 region of the 16S rRNA gene (Caporaso et al., 2011(Caporaso et al., , 2012Apprill et al., 2015;Parada et al., 2016). Library prep was conducted using iTru adapters for sequencing and a one-step PCR protocol with indexed primers, which is the current standard protocol used by the EMP (Caporaso et al., 2011(Caporaso et al., , 2012Apprill et al., 2015;Parada et al., 2016;Walters et al., 2016). Samples were sequenced on an Illumina MiSeq at the UAF Genomics Core Lab. Raw sequences were de-multiplexed using the Mr. Demuxy package (Cock et al., 2009). Demultiplexed sequences were run with mothur v1.40.0 on a high performance-computing cluster through UAF Research Computing Systems using a modified MiSeq standard operating procedure (Schloss et al., 2009). Operational taxonomic units (OTUs) were clustered at 100% similarity using the OptiClust option in mothur, taxonomy was assigned to OTUs using the SILVA 132 mothur formatted reference database with a bootstrap cutoff of 100%, and the samples in the resulting OTU table were normalized to 30,000 sequences and converted to relative abundances (Wang et al., 2007;Glöckner et al., 2017;Westcott and Schloss, 2017;Edgar, 2018). This OTU table (relative abundance of sequence reads per OTU for each sample) was used to conduct analyses of diversity and community structure of prokaryotes and to assess correlations with environmental variables.

Data Analyses
All statistical analyses were conducted using R, primarily using the vegan package (Oksanen, 2015;R Core Team, 2017). Diversity was quantified using the inverse Simpson index (1/λ) which reflects evenness and richness (Morris et al., 2014). Community structure was investigated via hierarchical clustering analysis with the Ward method (ward.D2) based on a Bray-Curtis dissimilarity derived from the OTU table. Cluster tests were used to assess the validity (heterogeneity, significance) and characteristics (silhouette, stability) of the hierarchical clusters. The heterogeneity of and significance between hierarchical clusters were investigated using Bray Curtis distances with the betadisper and adonis functions, respectively . OTUs that distinguished each cluster were identified via indicator taxa analysis using the indicspecies package (De Cáceres and Legendre, 2009;De Cáceres et al., 2010). The multipatt function was used to identify taxa specific to a single cluster, or indicative of combinations of clusters, with a significance of ≤0.001 and a strength of ≥0.900 (i.e., >90% probability of occurrence within a given cluster or combination of clusters (De Cáceres and Legendre, 2009). Given the complexity of the prokaryotic species concept and the lack of taxonomic resolution for many taxa, we used this indicator analysis to identify OTUs potentially indicative of certain habitat features and/or assemblages. We assessed the distribution of these indicator taxa at varying taxonomic levels, from phylum to genus, to characterize representative assemblages for each cluster.
Differences in the median values of environmental parameters among clusters were evaluated using Kruskal-Wallis tests followed by pairwise Wilcoxon Rank Sum tests with a Benjamini and Yekutieli (BY) correction (Yekutieli and Benjamini, 2002). Relationships between prokaryotic community structure and environmental variables were modeled using Constrained Analysis of Principle Coordinates, i.e., the capscale function in the R vegan package (Oksanen, 2015;R Core Team, 2017). The best model for capscale analyses was identified based on variance inflation factors for each environmental parameter and a forward and backward stepwise model selection via permutation tests on adjusted R 2 and P-values using the ordistep function. The significance of the environmental correlates yielded via Constrained Analysis of Principle Coordinates was investigated in R using the adonis function; correlation strengths were calculated using the cor function . All bioinformatic and statistical scripts can be found here on GitHub 1 .
In order to put Beaufort Sea surface sediment microbes into a broader geographical context, we qualitatively compared the 10 most abundant taxa at class level with those reported for global benthic surface sediments, global deep-sea sediments, and Norwegian-Arctic sediments (Zinger et al., 2011;Bienhold et al., 2016). We refer to the Arctic microbiome published by Bienhold et al. (2016) as the Norwegian Arctic microbiome because all but one dataset included in that study was generated from the Norwegian Arctic (Bienhold et al., 2016).

Prokaryotic Community Structure and Diversity
The total number of OTUs yielded was 295,115. No OTUs were removed (e.g., singletons or doubletons) because removals resulted in substantial changes in certain taxonomic groups (Supplementary Figure 1). Hierarchical clustering analysis identified four clusters, i.e., assemblages of prokaryotes which have significantly different (adonis: F = 12.6, P < 0.001) taxonomic composition (Figure 1). However, the silhouette test yielded an average cluster width of 0.11, indicating a high degree of similarity between clusters which suggests that the four assemblages may exist more so as a gradient rather than as discrete communities in particular locations. The number of samples within each cluster, herein referred to as assemblages A, B, C, and D, was as follows: A = 23, B = 35, C = 11, and D = 16. Assemblage A is most distinct, whereas C and D are most similar (Figure 1). The distribution of values for Inverse Simpson diversity differed for each assemblage, with the highest prokaryotic diversity exhibited in assemblage A, then B, D, and C.
Each assemblage is represented by sampling locations distributed across the study area (Figure 1), such that assemblages were not restricted to particular geographic areas that might reflect distinct environmental characteristics. However, when delving into the substructure of individual assemblages, patterns emerged at smaller spatial scales (Figure 2). Assemblage A contained two clusters, with one containing mostly Alaskan and Canadian Beaufort (AKCA) samples and some shallower (<350 m) Amundsen Gulf and Banks Island (AGBI) samples, and the other containing mostly AGBI samples with a few deeper (>350 m) AKCA samples. Assemblage B was clearly divided into AKCA and AGBI samples. Sub-structure for assemblage C was delineated by the Mackenzie River into the AK Beaufort shelf (west of the river) and the CA Beaufort with AG samples (east of the 1 https://github.com/amwalker/Beaufort-benthos river). Assemblage D divides into two groups with no obvious geographical demarcations.

Taxonomic Composition and Indicator Taxa
There were 350 prokaryotic families represented in this dataset; we focused on the 25 most abundant families (or in some cases, orders) to assess broad-scale patterns in community structure and compare relative abundance of these dominant taxa among individual assemblages. Abundant taxa consisting of multiple unclassified OTUs within a class or phylum, specifically Alphaand Gammaproteobacteria, were ignored. The 25 most abundant taxa accounted for 58% of the total sequence reads, and are listed in order of decreasing abundance in Figure 3. Familylevel composition differed among assemblages, with inverse relationships in abundance of particular taxa in assemblages A and C. Families that were relatively abundant in assemblage A also tended to be abundant in assemblage B, whereas assemblages C and D also showed similar trends. The families Nitrosopumilaceae, Woeseiaceae, NB1-j, Actinomarinales, Kiloniellaceeae, Subgroup 22, OM190, and AT-s2-59 were particularly abundant in assemblage A. Haliaceae, BD7-8, and Methyloligellaceae were most prominent in assemblage B. In assemblage C, Flavobacteraceae, Rhodobacteraceae, Rubritalaceae, and Psychromonadaceae were more abundant. Gammaproteobacteria Incertae sedis, Desulfobulbaceae, Sva1033, Anaerolineaceae, and Desulfobacteraceae were abundant in assemblage D.
Indicator taxa analyses revealed that assemblage A had the highest number of indicator OTUs (294), followed by C (61) Table 1). Taxa indicative of a combination of assemblages will be referred to as indicative taxa, to differentiate from assemblage-specific indicator taxa. Mirroring trends shown in Figure 3, assemblages A + B and C + D shared relatively high numbers of indicative OTUs. No indicative taxa were shared between assemblages A + C alone, but some were shared between A + B + C, suggesting B represents a transitional community with representatives from both of the distinct A and C assemblages. Assemblages B + D also shared very few indicative taxa, but more were found when combined with either assemblage A or C (A + B + D or B + C + D).
In light of the numerous uncultured and indicative OTUs, we limit our discussion to selected assemblage-specific indicator taxa, identified to the highest taxonomic resolution possible, for which published information could be used with confidence to infer function or distribution of a specific assemblage, and/or metabolic type or relationship to environmental variables (Figure 4). Indicator OTUs for assemblage A belonged to the class Acidobacteria (Subgroups 6,9,10,21,22,and 26), clades SAR202 and 324, and the families Thermoanaerobaculaceae, Magnetospiraceae (Magnetospira), Nitrosopumilaceae (Candidatus Nitrosopumilus), Nitrospiraceae (Nitrospira), and Scalinduaceae (Scalindua). Representative taxa under Thermoanaerobaculaceae belong to Acidobacteria subgroup 10 which was recently assigned to this family FIGURE 2 | Substructure of prokaryotic assemblages. This dendrogram highlights the substructure within each assemblage and the differing east to west gradients reflected by the substructure. Assemblage A was divided by location and depth, with one sub-cluster dominated by Amundsen Gulf (AG) samples combined with broader Beaufort (AKCA) slope samples deeper than 350 m. The other assemblage A sub-cluster was dominated by broader Beaufort slope samples with AG samples shallower than 350 m. Assemblage B divided cleanly between broader Beaufort and Amundsen Gulf samples. Assemblage C divided between the west of the Mackenzie River, AK samples, and east of the Mackenzie River, CA Beaufort and AG samples. Assemblage D divided into combinations of AG samples with either AK or CA samples with no obvious demarcations. (Dedysh and Yilmaz, 2018). Indicator taxa for assemblage C included representatives from the families Flavobacteriaceae (Polaribacter and Ulvibacter), Rhodobactereacae (Loktanella and Octadecabacter), Thiovulaceae (Cocleimonas), Thiotrichaceae (Sulfurimonas), and Desulfobulbaceae (Desulfobulbus) (Figure 4). Assemblage D indicator taxa were within the families Anaerolineaceae, Desulfobacteraceae, Desulfobulbaceae, orders PHOS-HE36 and SBR1031, and phyla Bathyarchaeota, and Schekmanbacteria. The only indicator OTU identified down to genus level here was Desulfococcus (Desulfobacteraceae), therefore other abundant genus-level OTUs that belonged to the families Desulfobacteraceae and Desulfobulbaceae, were also queried and further explored with the indicator taxa (Figure 4).

Environmental Correlates of Community Structure
Overall, the suite of environmental variables measured here explained relatively little (18%; Adj. R 2 = 0.18) of the The plot is arranged such that the most abundant taxa starts from the top and decreases toward the bottom. Note that Gammaproteobacteria Inc. sed., NB1-j, Actinomarinales, Sva1033, OM190, BD7-8, AT-s2-59, and UBA10353 are not classified to family, but order level. variation in prokaryote community structure among sites ( Figure 6). Nonetheless, results of the CAP analysis did suggest significant correlations (adonis: P < 0.001) with different habitat characteristics for each assemblage. Assemblage A was found solely on the continental slope in 23 samples ranging from 171 to 1,200 m water depth (median depth 350 m), thus generally consisting of "deep-sea" locations further offshore and below the photic zone. Sediments at these locations exhibited significantly higher δ 15 N and porosity, and lower values of Chl-a concentration and C:N; bottom waters were warmer and more saline (Figure 5). Assemblage A largely differentiated from the other assemblages on the CAP1 axis (Figures 6A,C). CAP1 was most strongly and positively correlated with depth, δ 15 N, and porosity and exhibited a moderate to weak negative correlation with Chl-a and TOC. Variation within assemblage A was most apparent on the CAP4 axis suggesting Amundsen Gulf and Banks Island samples (+CAP4) are more influenced by degraded and marine derived OM (higher δ 15 N and δ 13 C) whereas Beaufort slope samples (-CAP4) had higher porosity sediments with comparatively more phaeopigment and terrestrial OM ( Figure 6C).
Assemblage B consisted of 35 samples ranging between 36 and 350 m water depth (median depth 75 m), with a little more than half (19) located on the shelf and the remainder (16) on the slope. Median values of all environmental parameters for assemblage B were not significantly different than those for any of the other assemblages ( Figure 5). CAP analysis showed a separation of assemblage B largely along the CAP2, axis suggesting that stations in this assemblage are characterized by more degraded OM (positive correlation with phaeopigment and δ 15 N) and less influence from terrestrial organic matter (higher δ 13 C suggesting marine or ice-algal source; Dunton et al., 2006 ; Figures 5, 6A). Differences among samples within assemblage B were most apparent on both the CAP3 and CAP4 axes, indicating that Amundsen Gulf and Banks Island samples (-CAP3) were generally characterized by higher TOC and lower porosity, whereas samples from the FIGURE 4 | Proportion of assemblage-specific indicator taxa across prokaryotic assemblages. Indicator taxa are identified to the lowest taxonomic level possible. The plot is colored by functional groups which provide insight into characteristics associated with specific assemblages. *Indicates those taxa that were not yielded as indicators, but represent OTUs belonging to the same family as unidentified indicator taxa.
open shelf and slope (+CAP3) were characterized by higher porosity sediments with comparatively less TOC and more degraded phytodetritus (higher phaeopigment). Easternmost samples from the Amundsen Gulf (+CAP4) were separated from all other samples (-CAP4), with potentially sea-ice derived OM (higher δ 13 C), compared to more degraded marine and terrestrial OM (lower δ 13 C and higher phaeopigment) in the western Beaufort. Depth was also moderately or highly correlated with the CAP 3 and 4 axes, although this may reflect effects of other environmental variables that covary with depth.
The 11 samples represented by assemblage C were located solely in shallow shelf sediments ranging from 17 to 42 m water depth (median depth 21 m) with significantly lower values of δ 15 N, δ 13 C, and salinity. Assemblages C, B, and D were largely separated along the CAP2 axis, with assemblage C samples more correlated with fresh phytodetritus (higher Chl-a, lower δ 15 N) and terrestrial OM (lower δ 13 C). Samples within assemblage C varied along CAP3 axis which is most strongly correlated with TOC and porosity ( Figure 6B). Nearshore samples from the Alaskan Beaufort shelf (+CAP3) were characterized by lower TOC and higher porosity, whereas those from the Canadian Beaufort and Amundsen Gulf (-CAP3) had higher TOC and lower porosity.
Assemblage D consisted of 16 samples which were collected almost exclusively at shelf locations (=100 m; median depth 72 m), except for 3 deep sites (=200 m) offshore of the Colville River plume where high concentrations of Chl-a were found. Environmental characterization for and correlations with assemblage D were not well captured with parameters measured here (Figures 5, 6). Variation within assemblage D was most apparent on the CAP3 axis, which suggests that several samples had comparatively higher organic matter input and lower porosity.

Broader Geographical Context
Comparisons of the top 10 class-level taxa between Beaufort Sea, global benthic, deep-sea, and Norwegian-Arctic surface sediments indicate that all four sediment microbiomes share five major taxa: Gammaproteobacteria, Alphaproteobacteria, Deltaproteobacteria, Bacteroidia, and Actinobacteria (Figure 7). The Beaufort Sea sediment microbiome included taxa that were also highly abundant in either the combination of global benthic and deep-sea microbiomes (Plantomycetia) or the Norwegian-Arctic microbiome (Phycisphaera). Interestingly, four classes prevalent in the other three microbiomes were not among the 10 most abundant in Beaufort Sea sediments: Clostridia, Betaproteobacteria, Bacilli, and Acidobacteria. Nonetheless, Acidobacteria was FIGURE 5 | Distribution of environmental parameters for each assemblage. Each violin plot shows the distribution, within each assemblage, of the environmental parameters indicated as significant correlates in the ideal model for CAP analysis. Note that salinity is also included in this figure, though it was not run in CAP analysis due to the covariance inflation factor with depth. The symbols ✣ and indicate significant differences between un-matching symbols with other assemblages. δ 13 C values for ice, marine, and terrestrial OM sources was derived from Dunton et al. (2006). quite abundant in the Beaufort Sea, falling within the top 13 most abundant classes. Class Betaproteobacteria was not detected in Beaufort Sea surface sediments, whereas Clostridia and Bacilli were present but in relatively low abundances. Nitrososphearia, Verrucomicrobia, and Anaeronlinae were among the top 10 taxa of the Beaufort Sea sediment microbiome, but were not abundant in the other datasets.

DISCUSSION
We identified four distinct assemblages, A, B, C, and D, across a broad area in southern Beaufort Sea surface sediments. Assemblage A occurred only at sites on the upper continental slope with lower concentrations of more degraded OM and higher porosity sediments. A distinct community (assemblage C) was found in shallow shelf sediments, with higher OM content and terrestrial OM influence. The wide distribution of assemblage B across shelf and slope depths and overlap in taxonomic composition with assemblages A and C, may reflect a transition zone between the shelf and slope sites containing generalists found in both bathymetric zones. Assemblage D was found primarily at shelf depths, but included a few deeper slope sites in the westernmost part of the study area with high OM content. Indicator taxa with known metabolic functions and/or ecological context suggest that assemblage D was not merely comprised of generalists. In order to provide ecological context and better understand what distinguishes these assemblages, we examined taxa that were assemblage-specific and/or exhibited higher abundances associated to an individual assemblage together with environmental correlates.
Slope assemblage A may also reflect some of the unique biogeochemical characteristics on the Beaufort Sea slope. In the Arctic Halocline waters over the Canada Basin, SAR 202 have been implicated in terrestrial OM degradation, consistent with the relatively strong terrestrial signal exhibited throughout FIGURE 6 | Environmental parameters influencing prokaryotic assemblage structure. The ordination plots above depict the results of Constrained Analysis of Principal Coordinates. All included environmental parameters are significant correlates (P < 0.001) to the first four CAP axes (A-C) and, in total, explained 18% of the variation in prokaryotic community structure (R 2 = 0.18). The corresponding correlations are reported in the table for all axes.
the Southern Beaufort, even in the upper slope sediments (Figure 5; Colatriano et al., 2018). SAR202 has also been linked with deep-sea mud volcanoes and SAR324 with methane-rich hydrothermal plumes (Sheik et al., 2014;Coelho et al., 2016). Acidobacteria subgroup 10 was recently placed in the family Themoanaerobaculaceae which encompasses thermophilic obligate anaerobes isolated from marine hydrothermal vents and freshwater hot springs (Losey et al., 2013;Dedysh and Yilmaz, 2018). Several mud volcanoes on the Beaufort Sea slope release subsurface methane and produce temperatures ≥23 • C within a few centimeters of the sediment-water interface, creating suitable habitat for thermophilic anaerobes such as those within subgroup 10 (Paull et al., 2015;Gamboa et al., 2017;Lee et al., 2018;Paull, pers. comm.). Subgroup 10 is widespread across the study area, yet to date, mud volcanoes have only been reported on the Canadian Beaufort Slope. These microbes may be spore-forming, similar to other thermophilic anaerobes reported from Arctic marine sediments that produce spores that lie dormant until favorable conditions arise (Hubert et al., 2010;Müller et al., 2014). However, previous reports of thermo-anaerobic bacterial spores belonged to the phylum Firmicutes, and not Acidobacteria, and none of the current taxa within Acidobacteria subgroup 10 are known to be spore-forming (Hubert et al., 2010;Losey et al., 2013;Müller et al., 2014). Given their occurrence in Arctic sediments from other regions, it is more likely that bacteria in subgroup 10 inhabit a wider temperature range than previously recognized (Hoffmann et al., 2017).
Beaufort sediment mineralogical properties may be influenced by the abundant indicator OTUs within Magnetospiraceae, a family composed of magnetotactic bacteria. Magnetotactic FIGURE 7 | Top 10 abundant prokaryotes at class level exhibited in marine sediments. The colored grid above is arranged with surface sediment microbiomes on the y-axis, Beaufort, global benthic, Norwegian Arctic, and deep-sea (Zinger et al., 2011;Bienhold et al., 2016). The x-axis represents the rank from the most abundant taxa on the left to the least abundant taxa on the right exhibited by these microbiomes. The colors represent the class-level taxonomy of these prokaryotes.
bacteria are unique in that they biosynthesize a mineral nanocrystal which they can use as a "compass" to navigate toward favorable redox conditions found at the oxic-anoxic interface in marine sediments (Lefevre and Bazylinski, 2013). Though indicator OTUs were all uncultured, Magnetospira is the most abundant taxon identified to genus-level within this family. Members of the genus Magnetospira produce magnetite crystals, which are the predominant magnetic grain found in Canadian Beaufort sediments and are more prevalent on the slope Lefevre and Bazylinski, 2013;Gamboa et al., 2017).
Genera tightly linked to nitrification and anaerobic ammonia oxidation were relatively more abundant in slope assemblage A. Nitrospira are exclusively chemoautotrophic nitrite-oxidizing bacteria, oxidizing nitrite to nitrate, but have recently been found to perform complete nitrification, oxidizing ammonia to nitrate (Daims, 2014;Daims et al., 2015). Candidatus Nitrosopumilus encompasses ammonia-oxidizing archaea (AOA), which are aerobic nitrifiers that exhibit an extremely high affinity for ammonia, and can oxidize ammonia to nitrate, nitrite, or nitrous oxide (Stahl and de la Torre, 2012; Könneke et al., 2014;Martens-Habbena et al., 2015;Qin et al., 2016). AOA can perform mixotrophy and heterotrophy, but are largely chemoautotrophs (Stahl and de la Torre, 2012;Offre et al., 2013;Könneke et al., 2014;Qin et al., 2014). They also exhibit photoinhibition and increased autotrophic activity, which is consistent with the distribution of this assemblage at deeper sampling stations (Alonso-Sáez et al., 2012;Qin et al., 2014). In the pelagic zone of the Beaufort Sea, AOA are relatively more abundant in bottom waters on the slope and may be capable of urea degradation to fuel nitrification, as has been observed for Candidatus Nitrosopumilus sediminis in Arctic sediments near Svalbard (Alonso-Sáez et al., 2012;Park et al., 2012;Damashek et al., 2017). In the portion of our study area characterized by assemblage A, a peak in urea concentration occurs at 250-300 m water depth due to zooplankton and fish excretion, which may fuel the activity of AOA in this region (Simpson et al., 2008).
The genus Candidatus Scalindua encompasses anaerobic ammonia oxidizing (anammox) bacteria which convert ammonia to N 2 via autotrophic metabolism and can utilize nitrite, nitrate, metal oxides, and even oligopeptides and small organic molecules as electron acceptors (Penton et al., 2006;Van de Vossenberg et al., 2013). Anaerobic ammonia oxidation may be a major sink for nitrate in the marine system, and has been reported in Arctic shelf slope sediments in the Chukchi Sea, near Greenland, and in deep-sea sediments from the Arctic mid-ocean ridge (Rysgaard et al., 2004;Penton et al., 2006;Jorgensen et al., 2012). The predominance of anammox bacteria at the deeper slope sites occupied by assemblage A supports previous claims that anammox is more likely to occur in areas with less organic matter loading (Thamdrup and Dalsgaard, 2002;Rysgaard et al., 2004;Burgin and Hamilton, 2007;McTigue et al., 2016).

Generalist Assemblage (B)
While no indicator taxa were identified for assemblage B, one abundant taxon included members of the bacterial family Methyloligellaceae, dominated in this study by the genus Methyloceanibacter which contains known methylotrophs (Takeuchi et al., 2014;Vekeman et al., 2016a). Some species of Methyloceanibacter oxidize methanol, while others exhibit methane oxidation and include the methane monooxygenase gene in their genome (Takeuchi et al., 2014;Vekeman et al., 2016a,b). We have detected a relatively high abundance of the methane monooxygenase gene in these sediments through metagenomics analyses, although further analyses are needed to attribute these genes to a specific taxonomic group (Walker et al. unpublished data). Nonetheless, the presence of Methyloceanibacter coupled with the abundant distribution of the methane monooxygenase gene in methane-permeated Beaufort Sea sediments suggests that aerobic methane oxidation may be an active process in this region, although further research confirming this activity would be required. As there are multiple sources of methane in these sediments, methane oxidation in the Beaufort Sea is likely an important process constraining methane concentrations reaching the sea surface (Lorenson et al., 2016). However, no evidence for methanotrophy has yet been reported in the water column or subsurface sediments, despite chemical evidence that oxidation of methane derived from ancient sources may be occurring near the sediment-water interface (Hamdan et al., 2013;Sparrow et al., 2018).

Shallow Shelf Assemblage (C)
Indicator taxa for assemblage C were characteristic of a shallowshelf environment influenced by fresh phytodetritus, riverine input, and sea ice, in keeping with elevated Chl-a and significantly lower values of δ 13 C and salinity observed at assemblage C locations. Many of these taxa may be deposited to the seafloor in association with sinking particles. Flavobacteriaceae occur in conjunction with phytoplankton blooms and subsequent deposition of phytodetritus (Junge et al., 2002;Bienhold et al., 2012;Teeling et al., 2012;Hoffmann et al., 2017). Genera within this family are typically aerobic, though there are some facultative anaerobic heterotrophs capable of degrading complex polysaccharides (Teeling et al., 2012;McBride, 2014;Xing et al., 2015). Two of the genus-level indicator taxa found here, Ulvibacter and Polaribacter, are often particle attached and tightly connected with phytoplankton blooms and sea ice (Junge et al., 2002;Bowman and McCuaig, 2003;Brinkmeyer et al., 2003;Teeling et al., 2012;Xing et al., 2015). Ulvibacter is one of the first and most abundant genera to respond to a phytoplankton bloom following algal lysis, whereas Polaribacter is most abundant in later stages of an algal bloom and can use sulfatases to access highly sulfated algal material (Teeling et al., 2012(Teeling et al., , 2016Xing et al., 2015;Helbert, 2017). Polaribacter species generally include psychrophilic bacteria that are most commonly found in polar and temperate regions. Until recently, they had only been identified in sea ice and seawater, but are now described from Arctic sediments offshore of Svalbard (Gosink et al., 1998;Bowman et al., 2012;Li et al., 2014;Rapp et al., 2018).
A few of the indicator taxa for assemblage C suggest that light availability may be a major factor influencing community composition. Genomic analyses of Polaribacter species, along with Loktanella and Octadecabacter, indicator taxa from the Rhodobacteraceae family, contain various rhodopsins (lightdriven proton pumps) and other light responsive genes, indicating an affinity for well-lit shallower waters (Gonzalez et al., 2008;Boeuf et al., 2013;Vollmers et al., 2013;Wu et al., 2014;Xing et al., 2015;Dogs et al., 2017). Loktanella, described as anoxygenic photoheterotrophs (AAP), have been found in Beaufort Sea surface waters including the Mackenzie River outflow where they were positively correlated with Chla (Boeuf et al., 2013). AAP photosynthesize but also require oxygen for breaking down organic matter with light-harvested energy, and are most abundant in coastal waters and eutrophic estuaries (Waidner and Kirchman, 2007;Boeuf et al., 2013). Photoheterotrophy could be an advantage in such environments that are heavily influenced by terrestrial material, providing an energy source to fuel more costly metabolic processes such as the degradation of humic substances (Boeuf et al., 2013). Octadecabacter species are psychrophilic, aerobic heterotrophs with gas vesicles and are abundant in sea ice, though a few taxa have been isolated from seawater and marine sediments (Gosink et al., 1997;Brinkmeyer et al., 2003;Vollmers et al., 2013;Lee et al., 2014;Billerbeck et al., 2015). Octadecabacter spp. have also been linked to phytoplankton blooms and detrital aggregates in the Arctic and elsewhere (Grossart et al., 2005;Rapp et al., 2018;Koedooder et al., 2019).
The remaining indicator taxa for assemblage C belong to the bacterial families Thiotrichaceae, Thiovulaceae, and Desulfobulbaceae, all of which are involved in the marine sulfur cycle. Members of the Desulfobulbaceae family, which are anaerobic sulfate reducers, were also found in assemblage D and discussed further below (Kuever, 2014b). Thiotrichaceae and Thiovulaceae, represented here by Cocleimonas and Sulfurimonas, respectively, are sulfur oxidizers (Tanaka et al., 2011;Garrity et al., 2015;Han and Perner, 2015). Sulfur oxidizing bacteria are generally limited by the availability of necessary electron donors (typically reduced sulfur compounds) and acceptors (oxygen, nitrate) favored in anoxic and oxic conditions, respectively (Wasmund et al., 2017;Jørgensen et al., 2019). Adaptations exhibited by Sulfurimonas and Cocleimonas provide a competitive advantage in the oxidation of elemental sulfur and other sulfur intermediates (Tanaka et al., 2011;Pjevac et al., 2014;Han and Perner, 2015;Wasmund et al., 2017).

Anoxic Assemblage (D)
The indicator taxa for assemblage D suggest a greater extent of anoxic sediment within the upper 1 cm layer. All indicator OTUs for Assemblage D are strict anaerobes within the families Desulfobulbaceae, Desulfobacteraceae, Anaerolineaceae, PHOS-HE36, the order SBR1033, and the phyla Bathyarcheaota and Schekmanbacteria (Kuever, 2014a,b;Evans et al., 2015;Anantharaman et al., 2018;Yamada and Sekiguchi, 2018). Members of the archaeal phylum Bathyarcheaota and candidate bacterial phylum Schekmanbacteria consisted of uncultured OTUs with less than 85% matches to published sequences in BLAST. Though inferences based on this level of taxonomy are highly speculative, representatives from both phyla have been most commonly sequenced from subsurface marine sediments (Inagaki et al., 2003;Biddle et al., 2006;Sørensen and Teske, 2006;Evans et al., 2015;Anantharaman et al., 2018;Lopez-Fernandez et al., 2018;Winkel et al., 2018).
Anaerolineaceae are comprised of anaerobic heterotrophs with a fermentative metabolism and have been sequenced in anoxic, organic-rich habitats such as sludge digesters, hydrothermal vents, and glacial fjord sediments off the coast of Svalbard (McIlroy et al., 2017;Zeng et al., 2017;Yamada and Sekiguchi, 2018). They are typically abundant in biofilms and implicated as hydrogenotrophs and alkane degraders, though few studies have examined marine representatives (Siegert et al., 2011;Sinkko et al., 2013;Al-Thani et al., 2014;Liang et al., 2015;van der Waals et al., 2017). Similarly, bacteria in the order SBR1031 have been primarily investigated in waste-water treatment plants where they produce biofilms (Xia et al., 2016;Wang X. et al., 2018). Marine representatives have been found in oil-contaminated estuarine sediments, and increase in abundance in seasonal microbial mats in the North Sea (Obi et al., 2017;Cardoso et al., 2019). Bacteria from the family PHOS-HE36 have also been predominately studied in the context of sludge digesters and biofilms, often sequenced with bacteria in the Anaerolineaceae family, but have also been found in marine sediments near methane seeps, hydrothermal vents, and within microbial mats (Dabert et al., 2001;Fujii et al., 2002;Trembath-Reichert et al., 2016;Wang L. et al., 2018;Zinke et al., 2018;Yang et al., 2019).
Indicator taxa from assemblage D are also connected to the sulfur cycle, specifically sulfate reduction, which is the dominant process responsible for anaerobic degradation of organic matter in marine sediments (Wasmund et al., 2017;Jørgensen et al., 2019). The abundance of several taxa within the families Desulfobacteraceae and Desulfobulbaceae suggest that sulfate reduction is a prevalent in both shallow shelf assemblage C and anoxic assemblage D. At family level, Desulfobacteraceae taxa were more abundant in assemblage D, whereas Desulfobulbaceae taxa were more evenly distributed between assemblages C and D, albeit with specific genera that were more abundant in shallow assemblage C. These families have shown responses to different substrates in Arctic sulfidic sediments, suggesting the taxonomic differences between assemblages C and D reflect differences in substrate availability among sites (Dyksma et al., 2018;Müller et al., 2018).
Desulfocapsa and Desulfobulbus (Desulfobulbaceae) were more abundant in shallow shelf assemblage C, and have often been sequenced in Arctic fjord sediments alongside sulfuroxidizing bacteria (SOB) such as Sulfurimonas and Cocleimonas; which were also identified as indicator taxa for assemblage C (Zeng et al., 2017;Trivedi et al., 2018;Buongiorno et al., 2019). This repeated co-occurrence in Arctic marine sediments suggests an affinity for similar environmental/redox conditions, and further exemplifies subtle differences between sulfur cycling in assemblages C and D.
The more abundant SRB taxa in assemblage D, particularly SEEP-SRB1 and SEEP-SRB4 clades, are most commonly sequenced from cold seeps (Schreiber et al., 2010;Kleindienst et al., 2012;Ruff et al., 2015). SEEP-SRB have been implicated in coupling sulfate reduction to the anaerobic oxidation of non-methane and methane hydrocarbons, however they remain largely uncultured (Petro et al., 2019). Desulfococcus spp., though commonly found in sulfidic sediments, are also associated with marine seeps, and may be major players in anaerobic hydrocarbon degradation in association with seeps Orphan et al., 2001;Kleindienst et al., 2012Kleindienst et al., , 2014.
The majority of samples in this assemblage were located in shallow shelf sediments, some of which were characterized by heavy OM loading and low porosity, conducive to oxygen depletion. Samples collected on the slope in assemblage D exhibited some of the highest concentrations of Chl-a measured here, and may have received a recent pulse of fresh phytodetritus leading to shallow anoxic conditions in sediments. Moreover, multiple indicator taxa for assemblage D are also associated with biofilms, microbial mats, seep environments, or subsurface sediments, suggesting the presence of these features at some sites. Microbial mats have been reported on Beaufort Sea sediments, mostly surrounding slow gas seeps (Paull et al., 2011). The inherently patchy distribution of benthic seeps and microbial mats coupled with the varying conditions that foster anoxic sediments may explain the lack of definable substructure exhibited by assemblage D and weak correlations with environmental parameters in this study.

Broader Geographical Context
Benthic prokaryotic community substructure on the Beaufort Sea shelf and upper slope broadly reflected longitudinal gradients in OM quality and quantity, as well as the influence of bathymetry. While bathymetric effects were particularly apparent in distinguishing slope assemblage A, shallow shelf assemblage C, and generalist assemblage B, redox conditions likely played a role in differentiating assemblages C and D. These findings highlight the importance of assessing prokaryotic communities on a broad spatial scale, given that community structure reflected environmental variation occurring over scales of meters to kilometers. This raises the question of how prokaryotic community structure in the Beaufort Sea compares to other regions within the Arctic and beyond.
At class-level, 60% of the top 10 taxa from Beaufort Sea sediments were shared with those from the global, deep-sea, and Norwegian Arctic sediment microbiomes (Figure 7). The classes Anaerolineae, Nitrososphearia, and Verrucomicrobiae were abundant and unique to the Beaufort Sea surface sediment microbiome. The bacterial class Anaerolineae and archaeal class Nitrososphearia are typified by the families Anaerolinaceae and Nitrosopumilaceae which, as mentioned above, primarily contain anaerobic bacteria associated with heavy OM loading as was prevalent in anoxic assemblage D, and ammoniaoxidizing archaea prevalent in slope assemblage A, respectively. Not surprisingly, Nitrososphearia was not among the most abundant taxa in the global, deep-sea, or Norwegian Arctic sediment microbiomes. Though common in marine surface sediments and typically more abundant than ammonia oxidizing bacteria, archaea have been overlooked in previous studies solely focused on bacteria (Park et al., 2010;Zinger et al., 2011;Bienhold et al., 2016). Their relatively high abundance here highlights the under-appreciated role of archaea in marine surface sediments. If archaea are excluded, the bacterial taxon Thermoanaerobaculia joins the top 10 most abundant taxa in our study. This group encompasses the anaerobic, thermophilic organisms largely represented by subgroup 10 an indicator taxon for slope assemblage A. This class was not among the most abundant taxa in sediments globally, in the deep sea, or in the Norwegian Arctic. Thermoanaerobaculia together with Anaerolineae suggest comparatively shallower anoxic conditions in Beaufort Sea surface sediments, spanning both slope and shallow shelf locations.
Betaproteobacteria, Gemmatimonadetes, Bacilli, and Clostridia were among the 10 most abundant taxa from global benthic, deep-sea, and Norwegian-Arctic microbiomes, yet were not prominent in Beaufort sediments. Betaproteobacteria, which were not present even at low abundances in Beaufort Sea sediments, were discovered to be potential laboratory contaminants in other studies (Kulakov et al., 2002;Grahn et al., 2003;Laurence et al., 2014;Salter et al., 2014). On a global scale Zinger et al. (2011) observed Bacilli and Clostridia in higher abundances in coastal sediments vs. open water or deep-sea sediments. Bacilli and Clostridia can be indicators of urban fecal contamination in coastal watersheds, lagoons, and nearshore seawater (Wu et al., 2010;Dubinsky et al., 2012;Dai et al., 2018). As the Beaufort Sea is extremely remote from urban life, it is not surprising that these taxa were less prevalent. Clostridia was included in the Norwegian-Arctic microbiome and Bacilli in the deep-sea microbiome, which are also presumably more pristine than typical coastal areas, suggesting that there are bacteria in these classes which are not necessarily correlated with urban waste (Bienhold et al., 2016).
The prokaryotic community structure in the Beaufort Sea differed somewhat from the Norwegian-Arctic, even at the relatively coarse taxonomic level of class, suggesting benthic microbes may reflect regional differences in hydrography, biogeochemistry, and bathymetry of Arctic shelf systems (Carmack and Wassmann, 2006). However, samples from the Norwegian-Arctic were largely collected at depths > 1,000 m, so these regional differences need to be further explored particularly across wider depth ranges.

Synthesis
The occurrence of strictly anaerobic taxa in all assemblages identified here indicates that Beaufort Sea surface sediments exhibit widespread anoxia within the 0-1 cm depth horizon, consistent with oxygen penetration studies conducted in other Arctic sediments (Kostka et al., 1999;Jørgensen et al., 2005). A larger proportion of anaerobic taxa occurred in shallow shelf (C) and anoxic (D) assemblages, indicating that, in general, oxygen depletion in nearshore sediments is comparatively more widespread than offshore, particularly in areas where organic matter loading and Chl-a were higher and porosity lower. Although this is a common trend when comparing shelf to slope sediments, it may be of particular importance in the Beaufort Sea when viewed through the lens of climate change (Orcutt et al., 2011). An overall net increase of river runoff and terrestrial organic matter has already been established for the Mackenzie River, and is predicted to continue over the coming years (Holmes et al., 2012;Doxaran et al., 2015). Thus, we might expect an increase in abundance of anaerobic taxa in shallow shelf assemblage C, with potentially more stations resembling or overlapping with anoxic assemblage D. There may also be increased taxonomic overlap between microbial communities associated with the Mackenzie River outflow and shallow shelf assemblage C. In addition to riverine input, shifts in primary production and sea-ice extent will likely be most apparent in the shallow shelf sediments, as the majority of indicator taxa in assemblage C have been strongly correlated to both phytoplankton blooms and sea-ice microbial communities.
An increase in OM input, whether terrestrial or marine derived, will likely translate to an increase in OM burial and subsequent methanogenesis in the already methane-permeated Beaufort Sea sediments . Characterizing the microbes involved in both methane production/degradation and associated processes in these sediments is thus of particular interest. The high relative abundance of potential methanotrophs exhibited in this study could, in addition to highlighting the presence of methane throughout sediments in the region, indicate an important constraint on methane efflux that requires further investigation.
This study also highlights the vast array of uncultured and unclassified microbes in Arctic marine sediments, particularly exemplified by slope assemblage A, which was the most divergent and diverse and contained the highest number of uncultured taxa across all assemblages. Anoxic assemblage D also exhibited a comparatively high abundance of uncultured microbes with representatives which are more rarely studied in the context of the marine system. The lack of knowledge regarding these taxa points to a need for more experimental and culture-based studies to elucidate their role in biogeochemical processes, that would provide insights into the links between microbial community structure and benthic ecosystem function. Although taxonomic classification does not directly equate to function, known microbes from slope assemblage A highlight how taxonomic composition may provide a framework for future studies. For instance, research investigating the role of nitrification and anammox in the Beaufort Sea would be most beneficial in slope sediments where microbes involved in aerobic/anaerobic ammonia oxidation and nitrite oxidation were abundant, and could reveal useful information about regional nitrogen cycling and upwelling dynamics. Additionally, microbial taxonomy can provide insights into environmental conditions that are not typically measured with surveys of larger benthic epi-and infaunal organisms, such as sediment oxygen penetration and the presence of biofilms, microbial mats, and/or seeps seemingly reflected by the taxonomic composition of anoxic assemblage D. Overall, these insights provided by prokaryotic community structure emphasize the value of incorporating microbial surveys more broadly into field sampling programs in the North American Arctic and elsewhere. This study establishes a baseline for Beaufort Sea benthic microbes, however without long-term datasets we cannot assess how these communities will respond to and reflect shifts related to climate change or other disturbances such as energy resource exploration.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA451041.

AUTHOR CONTRIBUTIONS
AW was involved in field work for sample collection, conducted all sample processing and data analysis, and took a lead role in the interpretation of results and manuscript preparation. ML contributed funding, advised on data analysis and interpretation, and contributed to manuscript preparation. SM conducted field work to collect samples and environmental data, provided funding, and contributed to project planning, interpretation of results, and manuscript preparation. All authors contributed to the article and approved the submitted version.

FUNDING
Sample sequencing, student stipend, and tuition were funded by the US Department of the Interior, Bureau of Ocean Energy Management (BOEM) through the Coastal Marine Institute graduate student fellowship program (Grant# M16AC00004). Stipend and tuition support for AW was provided in part by the North Pacific Research Board (Project #1303). Support for the use of the UAF Genomics Core lab was provided by the Institutional Development Award (IDeA) from the National Institute of General Medical Sciences of the National Institutes of Health under grant number 2P20GM103395. The content of this manuscript is solely the responsibility of the authors and does not necessarily reflect the official views of the NIH. This work was supported via in-kind match from the high-performance computing and data storage resources operated by the Research Computing Systems Group at the University of Alaska Fairbanks Geophysical Institute.

ACKNOWLEDGMENTS
Samples for this project were collected opportunistically through field work funded by BOEM and Department of Fisheries and Oceans Canada. Sincerest gratitude goes out to Dr. Holly Bik for providing the bioinformatics and genomics foundation needed for this project to happen. Appreciation also goes to Dr. Eric Collins for ongoing bioinformatic and statistical assistance. A huge thank you to Shannon MacPhee for coordinating and collaborating with us, without you we would not have had such an expansive dataset from the Canadian Beaufort Sea. The majority of environmental variables incorporated in this study were provided through the hard work of Julia Dissen.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2021.581124/full#supplementary-material Supplementary Figure 1 | Distribution of Methyloceanibacter at different levels of sequence removal. The pie chart above illustrates an example, using the genus Methyloceanibacter, of how relative abundances change with removal of singletons and doubletons compared to no removal (none). The difference in abundance is attributed to the removal of OTUs that would likely not have been removed if clustered at 97% similarity as they belong to the same genus. This is an important consideration when using amplicon sequence variants (ASVs), and software which automatically remove singletons and doubletons such as DADA2. Supplementary Table 1 | Indicator taxa for each prokaryotic assemblage. Taxa are arranged from most abundant starting at the top to least abundant toward the bottom. Those without a semi-colon indicate OTUs which could not be IDed past the taxonomic level provided.