16S and 18S rRNA Gene Metabarcoding Provide Congruent Information on the Responses of Sediment Communities to Eutrophication

Metabarcoding analyses of bacterial and eukaryotic communities have been proposed as efficient tools for environmental impact assessment. It has been unclear, however, to which extent these analyses can provide similar or differing information on the ecological status of the environment. Here, we used 16S and 18S rRNA gene metabarcoding to compare eutrophication-induced shifts in sediment bacterial and eukaryotic community structure in relation to a range of porewater, sediment and bottom-water geochemical variables, using data obtained from six stations near a former rainbow trout farm in the Archipelago Sea (Baltic Sea). Shifts in the structure of both community types were correlated with a shared set of variables, including porewater ammonium concentrations and the sediment depth-integrated oxygen consumption rate. Distance-based redundancy analyses showed that variables typically employed in impact assessments, such as bottom water nutrient concentrations, explained less of the variance in community structure than alternative variables (e.g., porewater NH4+ inventories and sediment depth-integrated O2 consumption rates) selected due to their low collinearity (up to 40 vs. 58% of the variance explained, respectively). In monitoring surveys where analyses of both bacterial and eukaryotic communities may be impossible, either 16S or 18S rRNA gene metabarcoding can serve as reliable indicators of wider ecological impacts of eutrophication.

Metabarcoding analyses of bacterial and eukaryotic communities have been proposed as efficient tools for environmental impact assessment. It has been unclear, however, to which extent these analyses can provide similar or differing information on the ecological status of the environment. Here, we used 16S and 18S rRNA gene metabarcoding to compare eutrophication-induced shifts in sediment bacterial and eukaryotic community structure in relation to a range of porewater, sediment and bottom-water geochemical variables, using data obtained from six stations near a former rainbow trout farm in the Archipelago Sea (Baltic Sea). Shifts in the structure of both community types were correlated with a shared set of variables, including porewater ammonium concentrations and the sediment depth-integrated oxygen consumption rate. Distance-based redundancy analyses showed that variables typically employed in impact assessments, such as bottom water nutrient concentrations, explained less of the variance in community structure than alternative variables (e.g., porewater NH 4

INTRODUCTION
Assessing the ecological integrity of benthic habitats is crucial to marine ecosystem management (Fernandes et al., 2001), with seafloor monitoring efforts having traditionally relied on morphological inventories of macrofauna (≥0.5 mm size fraction) and associated indices (Lejzerowicz et al., 2015;Pawlowski et al., 2016). While producing such inventories requires considerable time, taxonomic expertise and finances, environmental DNA (eDNA) metabarcoding has emerged as a high-throughput and low-cost option for characterizing the structure and composition of sediment communities . Comparisons of morphology-based and metabarcoding approaches for the characterization of sediment macrobenthos have demonstrated that the latter can function as a promising complement to traditional surveying methods (Aylagas et al., 2014(Aylagas et al., , 2016Cordier et al., 2017;Lobo et al., 2017;Cahill et al., 2018;Stoeck et al., 2018;Clark et al., 2020;Frontalini et al., 2020). By relying on total DNA extraction from samples, eDNA metabarcoding also enables investigations of microbial and meiofaunal (≤0.5 mm) communities that are extremely challenging to study by traditional means, providing a more comprehensive snapshot of ecosystem health than macrofaunal surveys alone (Chariton et al., 2015;Pawlowski et al., 2016;Stoeck et al., 2018;Frontalini et al., 2020).
Efforts to explore the use of metabarcoding in benthic ecosystem monitoring have involved analyses of both sediment eukaryotic (18S rRNA gene) and bacterial (16S rRNA gene) communities (reviewed in Pawlowski et al., 2018). In the context of aquaculture-induced organic enrichment, this work has focused on foraminifera (Pawlowski et al., 2014(Pawlowski et al., , 2016He et al., 2019) and ciliates , as well as general analyses of eukaryotic (Chariton et al., 2015) and bacterial communities (Dowle et al., 2015;Fodelianakis et al., 2015;Stoeck et al., 2018;Verhoeven et al., 2018;Moncada et al., 2019). The use of metabarcoding as a monitoring tool has also been explored with reference to localized disturbances caused by dredging (Zhang et al., 2017) and off-shore drilling (Lanzén et al., 2016;Laroche et al., 2017;Frontalini et al., 2020). These studies have shown that bacterial and eukaryotic metabarcoding analyses can serve as sensitive methods to detect anthropogenic impacts on sediment habitats, including both short-term responses and longterm effects on community resilience and stability. A recent cross-laboratory comparison of eDNA metabarcoding results produced using a standardized protocol has also demonstrated a high degree of reproducibility between individual metabarcoding data sets, which is essential to the successful deployment of eDNA-based monitoring methods (Dully et al., 2021).
Despite the promising results reported in previous studies and the increasing affordability of large-scale sequencing surveys, neither 16S nor 18S rRNA gene metabarcoding have yet been adopted as a routine component of monitoring programs (but see Lefrançois et al., 2018). While using both techniques in tandem would provide particularly comprehensive data on the status of sediment communities, this is likely to be unfeasible due to monitoring programs being subject to strict time and financial limitations (Borja and Elliott, 2013;Borja et al., 2016). It also remains unclear to what extent these methods can offer similar or contrasting information on the status of seafloor habitats, particularly with reference to identifying physical and chemical variables that are likely to drive community shifts at multiple trophic levels. Although several studies have correlated water column and/or sediment geochemical variables with metabarcoding data (Pawlowski et al., 2014;Chariton et al., 2015;Dowle et al., 2015;Fodelianakis et al., 2015;Lanzén et al., 2016;Forster et al., 2018;He et al., 2019;Moncada et al., 2019), there is a lack of research on this topic involving cross-comparisons of prokaryotic and eukaryotic communities. In cases where both types of metabarcoding analyses have been employed, shifts in community structure have either not been correlated with geochemical data or only a limited set of measurements has been used (La Rosa et al., 2001;Zhang et al., 2017;Keeley et al., 2018;Stoeck et al., 2018). Addressing this absence of information, therefore, is essential to establishing a full understanding of how metabarcoding can be best employed as a tool for marine ecosystem monitoring.
Here we used a field transect approach to compare 16S and 18S rRNA gene metabarcoding as tools to obtain insights into the impacts of aquaculture-induced eutrophication on the structure and composition of sediment communities in the coastal Archipelago Sea (Baltic Sea, Finland). In addition, we determined whether shifts in the structure of these communities were correlated with a shared or divergent set of environmental variables, and whether variables routinely included in benthic monitoring programs could reliably predict community responses to eutrophication. The data show that eutrophication-associated shifts in bacterial and eukaryote community structure are to a large extent linked to a common set of variables, suggesting that community changes detected by 16S or 18S rRNA gene metabarcoding are likely to reflect wider ecological responses to increased organic loading. Our findings additionally demonstrate that the sensitivity of monitoring approaches could be improved through more carefully designed protocols for the collection of environmental metadata, especially with reference to understanding the impacts of eutrophication at multiple levels of biological organization.

Study Location and Sampling Strategy
Six stations (S1-S6) were sampled in the vicinity of a former rainbow trout farm adjacent to Haverö Island, Archipelago Sea, Baltic Sea (Figure 1 and Supplementary Table 1). The farm was operational between 1987 and 2008 (Jokinen et al., 2018), after which it has been intermittently used for storing live fish. Sites S1-S3 were situated within a basin where the farm was located, S4 on a sill, S5 on the seaward side of the sill and S6 was included as a reference site outside the direct influence of the farm (Figure 1). At each site, water column variables were measured with a conductivity-temperature-depth (CTD) profiler (SBE 19 SEACAT, Sea-Bird Electronics Inc., Washington, United States), and sediment cores were collected with a GEMAX twin corer (internal ∅ 9 cm) in September 2017. Separate cores were collected for geochemical analyses (of porewater and bulk sediment), and for microbiological analyses.

Porewater Sampling and Analysis
Porewater samples for geochemical analyses were taken immediately upon retrieval of sediment cores at 2-cm depth intervals from 0 to 20 cm depth, using Rhizon samplers (Rhizosphere Research Products, Wageningen, Netherlands) (Seeberg-Elverfeldt et al., 2005). Two cores were taken for Rhizon samples, one for nutrient analysis and the other for analysis of hydrogen sulfide (H 2 S). From the first core, subsamples for analysis of dissolved nitrogen species were stored at -20 • C, while subsamples for inductively coupled plasma mass spectrometry (ICP-MS) analysis of phosphorus were immediately acidified with 1M nitric acid (HNO 3 ) and stored at 4 • C. Samples for H 2 S (from the second column) were prepared as described in Jilbert et al. (2018) with pre-addition of 1 ml of 2 M zinc acetate to the sampling syringe. A third core was subsampled with mini-cores (internal ∅ 2.3 cm; n = 1 for site S1 and n = 3 for S2-S6) for determination of the porewater oxygen profile.
Porewater nitrate (NO 3 − ), nitrite (NO 2 − ), and ammonium (NH 4 + ), concentrations were measured at Tvärminne Zoological Station, Hanko, Finland. NO 3 − and NO 2 − concentrations were measured using an autoanalyser equipped with a cadmium reduction column (Aquakem TM 250, Thermo Fisher Scientific TM , Waltham, MA, United States; detection limit of 0.2 µM), according to standard methods (Finnish Standards Association methods SFS 3030 and 3032). NH 4 + concentrations were measured manually with a detection limit of 0.14 µM, according to methods described in Koistinen et al. (2018). Total porewater phosphorus (P) concentrations were determined at Geo Lab, Utrecht University, the Netherlands, using a Thermo Scientific TM Element 2 TM ICP-MS (P determined as 31 P). Replicate analyses with reference to in-house standards indicated that the relative error for analyses of porewater nutrient concentrations was <5% in all cases.
Porewater H 2 S concentrations were determined by spectrophotometry as described in Jilbert et al. (2018). An acidic solution of ferric chloride (FeCl 3 ) and N,N-dimethyl-pphenylenediamine was added directly to sample vials. The zinc sulfide (ZnS) precipitate formed during sampling complexes S as methylene blue, allowing spectrophotometric analysis at 670 nm (Cline, 1969;Reese et al., 2011). H 2 S concentrations were calibrated against a series of standard solutions of sodium sulfide (Na 2 S·3H 2 O), whose S concentrations were determined by iodometric titration.
Inventories of porewater NH 4 + , dissolved P and H 2 S were calculated for the uppermost 20 cm of the sediment column, including bottom water concentrations. Depth-integrated concentrations are likely to provide more detailed information on biogeochemical processes taking place in the sediment than spot measurements of bottom water or surface sediment concentrations. The concentrations (µmol cm −2 ) are expressed as depth-integrated values (Inv x ): Where is the sediment porosity and Conc. (x) is the porewater concentration of NH 4 + , P, or H 2 S at each sediment depth interval or in bottom water (represented by a value of -1 in the above equation). The equation is modified from Canfield (1989).
Porewater oxygen profiles were measured with an oxygen microsensor (OX-100; Unisense A/S, Aarhus, Denmark). The microsensor was two-point calibrated in 100% air-saturated filtered seawater, collected from the study site, and in anoxic solution containing sodium ascorbate and NaOH (both at 0.1M). Measurements were carried out at depth intervals of 100 µm. Depth-integrated sediment pore water oxygen consumption rates were calculated from oxygen microprofile data using a least squares fitting routine in the PROFILE software package 1 (Berg et al., 1998).

DNA Isolation, PCR and Sequencing of 16S and 18S rRNA Genes
DNA was isolated from 0.25 g of surface sediment (top 1 cm; five biological replicates for each sampling site, each from a separate core) using a Powersoil R DNA isolation kit (MO BIO, Carlsbad, CA, United States) according to the manufacturer's instructions. A kit control (n = 1) with no added sediment was prepared. PCR was used to amplify 16S rRNA gene (V1-V3 regions) and 18S rRNA (V4 region) gene sequences ( Table 1). Each PCR contained 0.5 µl of template DNA, 1× Phusion TM Flash High-Fidelity PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, United States), and 0.5 µM of forward and reverse primer mixes, made up to a total volume of 25 µl with sterile ultrapure water. Technical duplicates and no-template controls were included in each PCR run. Cycling conditions were 98 • C for 10 s, followed by 20-22 or 28-32 cycles (16S and 18S rRNA gene, respectively) of 98 • C for 10 s and 72 • C for 15 s, and a final extension step at 72 • C for 1 min. PCR products were visualized following electrophoresis on 1% (w/v) agarose gels, after which the duplicates were pooled. Further sample processing was carried out by the DNA Sequencing and Genomics Laboratory, Institute of Biotechnology, University of Helsinki 2 . Custom barcodes for sample demultiplexing (Somervuo et al., 2018) were attached in a second PCR step and the obtained PCR products were purified using Agencourt R AMPure R XP magnetic beads (Beckman Coulter, CA, United States), quantified and pooled in iso-molecular quantities (Salava et al., 2017), followed by sequencing on the Illumina MiSeq platform. 1 | Primers and primer mixes used for the PCR amplification of 16S rRNA and 18S rRNA genes from DNA extracted from sediments collected from the Archipelago Sea (Baltic Sea).

Sequence Processing and Bioinformatics
Following the removal of MiSeq adapter and barcode sequences, further sequence processing was carried out using micca v. 1.6.2 3 (Albanese et al., 2015). Paired-end reads were merged using "micca mergepairs" ( Table 2; Rognes et al., 2016) and primer sequences were trimmed using "micca trim" (Martin, 2011). Quality filtering was performed using "micca filter" (filtering parameters in Table 2; Rognes et al., 2016). The command "micca otu" was used for chimera filtering and OTU clustering, with the clustering step employing a de novo greedy clustering algorithm with a 97% similarity threshold (parameters -d 0.97 -c) (Westcott and Schloss, 2015;Rognes et al., 2016). Taxonomic assignments were carried out using "micca classify" and RDP Classifier v. 2.11 (16S rRNA gene sequences; Wang et al., 2007) or the SILVA 132 database (18S rRNA gene sequences, majority taxonomy mapping file with seven taxonomic levels; Quast et al., 2013). Additional data processing was carried out using R v. 4.0.2 (R Core Team, 2020). To improve between-sample comparability, two replicates that markedly deviated from other replicates in terms of sequencing depth or observed OTU counts were omitted from the 18S rRNA gene sequence data set. Further replicates were randomly discarded to obtain a balanced design 3 https://compmetagen.github.io/micca/ for statistical analysis, with the final 18S rRNA gene sequence data set retaining four replicates for each site. For the 16S rRNA gene sequence data set, five replicates were retained for each site. Rarefaction curves for the resulting data sets are provided in Supplementary Figure 1. Unclassified phyla and 16S rRNA gene sequences annotated as Chloroplast (class level) or Mitochondria (family level) were removed. Within-station relative abundances of sequences unclassified at the domain level were <0.001% for the 16S rRNA gene sequence data set and <0.01% for the 18S rRNA gene sequence data set, respectively, with similar values observed for each sampling station. Removing unclassified, chloroplast and/or mitochondrial sequences had no major influence on sample clustering, as indicated by comparisons of non-metric multidimensional scaling (nMDS) ordinations  Figure 2). Prior to statistical analysis, both data sets were subjected to a denoising step using a 5% prevalence threshold (Callahan et al., 2016 ; Supplementary Figure 3). Prevalence is defined as the number of samples in which a taxon appears at least once (Callahan et al., 2016).
To explore correlations between community structure and sediment, sediment porewater and water column geochemical variables (Table 3), linear dependencies between variables were first identified by a principal component analysis (PCA) using centered and scaled data (Figure 2), and by inspecting variance inflation factors (VIFs). Two sets of up to five explanatory variables (Table 4) were then selected for distancebased redundancy analyses (db-RDA; Legendre and Anderson, 1999;Ramette, 2007) performed with vegan (Oksanen et al., 2020). Distance-based redundancy analyses were selected as the ordination method based on the inspection of primary axis lengths using detrended correspondence analyses (Lepš and Šmilauer, 2003). Variables in the first model (Model 1) were selected due to their minimal multicollinearity (VIFs of ≤6), while the second model (Model 2) was limited to variables similar to those often measured during environmental monitoring surveys (see e.g., the HELCOM Monitoring Manual 5 ) and had VIFs of ≤13.7 (Table 4). Both models were run using CLRtransformed 16S or 18S rRNA gene sequence data, resulting in a total of four model runs. Environmental variables were projected onto the ordination space using the envfit() function in vegan  (Oksanen et al., 2020), with fitting performed on linear scores of the ordination axes. Following global significance tests for the ordinations, the variance explained by each variable was compared using permutation tests with the remaining variables included as covariates, using the anova() function in vegan (999 permutations; Oksanen et al., 2020).

Bacterial Community Structure and Composition
Good's coverage estimates of 96.13-98.70% were obtained for the raw 16S rRNA gene sequence data set (see Supplementary Figure 1A for rarefaction curves). Following quality filtering and denoising, a total of 1,089,217 reads (clustered into 4,499 OTUs) were retained in the data set. Sediment bacterial community structure differed between sampling stations, as demonstrated by metabarcoding analysis of a total of 30 microbial communities (1-way global PERMANOVA using Aitchison distances: pseudo-F 5 , 24 = 6.559, p < 0.001) (Figure 3A). Post hoc pairwise comparisons showed that all six stations were distinct in terms of bacterial community structure (PERMANOVA with Benjamini-Hochberg correction: p = 0.014 for all tests). A comparison of inter-site variation in multivariate dispersions showed these differences to be due to a location effect rather than a dispersion effect (PERMDISP: pseudo-F 5 , 24 = 1.628, p = 0.203).
Comparing OTU relative abundances showed sediments from all sites to be dominated by 16S rRNA gene sequences from the phylum Proteobacteria, followed by Bacteroidetes (Figure 4A; for class-level relative abundances, see Supplementary Figure 4). The relative abundances of Proteobacteria were comparatively stable along the transect (S1-S5; values of 34-47%), with slightly lower values (30-32%) within the reference site (S6) ( Figure 4A). However, there were up to three-fold differences in the relative abundances of the four next most dominant phyla, including members of the Bacteroidetes, Acidobacteria, Chloroflexi, and Ignavibacteriae ( Figure 5A). For example, sediments within the site nearest to the fish farm (S1) showed a high relative abundance (>25%) of sequences from the phylum Bacteroidetes (classes Flavobacteriia and Bacteroidia; Supplementary Figure 4) with low relative abundances of Acidobacteria and Chloroflexi in comparison with sites farther from the farm (S4-S6) ( Figure 5A). The relative abundances of Ignavibacteriae were lower in S1 in comparison with other sites while in stations within the basin in the vicinity of the fish farm (i.e., S2 and S3), relative abundances of Acidobacteria were also lower than in S4-S6 ( Figure 5A). In contrast with this variation in relative abundances, no clear intersite variation in bacterial diversity was observed along the study transect (S1-S5) or when comparing sites to reference station S6, with mean Shannon's diversity ranging from 6.25 (S6) to 7.13 (S4) ( Table 5).

Eukaryotic Community Structure and Composition
Good's coverage estimates of 99.98-100% were obtained for the raw 18S rRNA gene sequence data set (see Supplementary Figure 1B for rarefaction curves). A total of 6,902,650 reads (clustered into 548 OTUs) were retained in the data set following quality filtering. Similar to the bacterial data, metabarcoding analysis of 24 sediment eukaryote communities demonstrated significant differences between stations (1-way global PERMANOVA: pseudo-F 5 , 18 = 2.415, p < 0.001; Figure 3B). Each station had a distinct eukaryote community structure (post hoc pairwise PERMANOVA with Benjamini-Hochberg correction: p = 0.036 for all tests). In contrast with the bacterial data, significant inter-site variation in multivariate dispersion was observed (PERMDISP: pseudo-F 5 , 18 = 9.260, p < 0.001). Sediment eukaryote communities at site S4 exhibited a higher degree of multivariate dispersion than communities at S1, S3, and S5, as shown by a post hoc Tukey's HSD test (Figure 6). However, no significant differences in multivariate dispersion (p > 0.05) were detected between site S6 and sites along the study transect (S1-S5; see Figure 6). Comparing OTU relative abundances showed sediments from all sampling stations to be dominated by unicellular eukaryotic taxa (Dinoflagellata followed by Protalveolata and Ochrophyta) (Figure 4B). While the relative abundances of Dinoflagellata appeared comparable between stations, nearly two-fold differences in overall eukaryotic diversity (Shannon's diversity) were observed between sites (site-specific means of 1.00-1.88; Table 5). For example, Metazoa ( Figure 5B) exhibited a significant decline in Shannon's diversity in sites near the fish farm (S1-S3;x = 1.17 ± 0.08) compared to sites farther from the farm (S4-S6;x = 1.89 ± 0.11) (two-sample Wilcoxon rank sum test, W = 11, p < 0.001). Metazoa detected in sediments from site S1 included copepods (Calanoida) and nematodes (Monhysterida), with these taxa also detected in other sites ( Figure 5B). Examples of taxa encountered in sites S4-S6 but either absent or present at low relative abundances (<0.2%) from S1-S3 included Enoplida (nematodes), Haplotaxida (annelids), Homalorhagida (kinorhynchs), and Macrostomida (flatworms) (Figure 5B).

Correlations Between Community Structure and Environmental Variables
Sediment depth-integrated O 2 consumption rates were two orders of magnitude higher in the basin (sites S1-S3; rates up to 0.016 µmol −1 cm −2 s −1 ) than within other sites (rates up to 0.006 µmol −1 cm −2 s −1 ) ( Table 3). Oxygen concentrations in the bottom water were also up to two times lower in sites near the farm than in other sites (minimum of 80.1 µmol L −1 in S1). Concentrations of porewater nutrients (NH 4 + , NO x , and P) and sediment organic carbon (C org ) were high in S1-S3, with the NH 4 + inventory in S1 (18.7 µmol cm −2 ) being approximately an order of magnitude higher than those in S4-S6. Differences in bottom water nutrient concentrations were also observed. Bottom water P concentrations were relatively low in sites S2-S6 (2.3-6.3 µmol L −1 ), with a high concentration detected at S1 (42.3 µmol L −1 ). Variation in bottom water NO x concentrations was less pronounced, with the highest concentration (14.3 µmol L −1 ) observed at S1. While the sediment within site S4 consisted of fine sand (median grain size of 93.6 µm), all other sites were muddy (median grain sizes of 3.1-13.1 µm). Further information on the environmental conditions within each sampling station is given in Table 3.
Correlations between site-specific environmental conditions and sediment bacterial or eukaryotic community structure were explored using distance-based redundancy analyses (db-RDA; see Table 6) performed using two separate sets of up to five environmental variables (Models 1 and 2; details on model construction in Statistical Analysis and Table 4). The first db-RDA model (Model 1) explained 58% (global permutation test: pseudo-F 5 , 24 = 6.559, p < 0.001) and 40% (global permutation test: pseudo-F 5 , 18 = 2.415, p < 0.001) of the variance in sediment bacterial and eukaryotic community structure, respectively (Figures 7A,B). As shown by permutation tests for individual db-RDA explanatory variables (with the remaining variables used as covariates), in Model 1 the observed variation in bacterial community structure was strongly correlated with porewater NH 4 + concentrations (pseudo-F = 7.230, p B−H = 0.005), followed by sediment depth-integrated O 2 consumption rates and distance to the farm ( Table 6). A near-significant correlation between bacterial community structure and the sediment C org :N tot ratio was observed (pseudo-F = 1.924, p B−H = 0.065), while the weakest correlation (pseudo-F = 1.685, p B−H = 0.09) corresponded to sediment grain size ( Table 6). In contrast with bacterial communities, all variables included in Model 1 were significantly correlated with variation in eukaryotic community Frontiers in Marine Science | www.frontiersin.org FIGURE 5 | Operational taxonomic unit relative abundances (%) corresponding to subsets of (A) bacterial and (B) eukaryotic communities within sediments collected from a total of six sampling sites in the Archipelago Sea (Baltic Sea, Finland). The bacterial community subset corresponds to the top four most abundant phyla, aside from members of the Proteobacteria. The eukaryotic community subset corresponds to metazoa. structure ( Table 6). The environmental variables most strongly correlated with variation in eukaryotic community structure corresponded to the sediment C org :N tot ratio, porewater NH 4 + concentrations and sediment depth-integrated O 2 consumption rates, followed by distance to the farm ( Table 6).
( Figures 7C,D). Significance values of all tested variables for both models and effect sizes (pseudo-F values) are shown in Table 6.
In contrast to C org :N tot (in Model 1), significant correlations (p B−H < 0.05) were observed between the sediment C org content and both bacterial and eukaryotic community structure. Shifts in the structure of both community types were also correlated with bottom water NO x , P, and O 2 concentrations. Effect sizes within the bacterial data set often exceeded those in the eukaryotic data set. For example, an over fourfold difference in effect sizes was observed in relation to sediment porewater NH 4 + inventories (see Table 6).

DISCUSSION
16S and 18S rRNA gene metabarcoding have been proposed as alternatives to morphology-based macrofaunal inventories FIGURE 6 | Distances to the group centroid for CLR-transformed eukaryotic OTU abundance data. Data are shown for a total of six sampling sites in the Archipelago Sea (Baltic Sea, Finland). Different letters above the boxes indicate significant differences between sampling sites (Tukey's HSD test, p < 0.05).
that have been traditionally used in benthic monitoring surveys (reviewed by Pawlowski et al., 2018). To date, little information has been available concerning the ability of these methods to provide similar versus contrasting information on ecosystem health. Here, we employed a field transect study to compare shifts in bacterial and eukaryotic community structure in relation to 18 sediment and water column variables along an aquaculture-induced eutrophication gradient near a former rainbow trout farm (active during 1987-2008Jokinen et al., 2018) in the Archipelago Sea (Baltic Sea, Finland). Correlations between community structure and environmental variables were explored using two separate models including environmental variables selected either on account of their low multicollinearity or because they are often measured during monitoring surveys.
A key feature of the data was that the site nearest to the farm (S1) showed high relative abundances of taxa belonging to the phylum Bacteroidetes in parallel with a reduction in the abundance of other bacterial taxa (members of the phyla Acidobacteria and Chloroflexi), compared with other stations. Relative abundances of Ignavibacteriae were also low in site S1 compared to other sites. Moreover, Acidobacteria in sediments from sites S2 and S3 exhibited low abundances in comparison with sites farther from the farm (S4-S6). These shifts in bacterial community structure are in agreement with findings reported by Quero et al. (2020), who conducted a 10-months study into the effects of on-going sea bass and sea bream farming on the structure of sediment bacterial communities in Southern Sicily. Similar to our study, the authors reported increased relative abundances of Bacteroidetes in fish farming-impacted sediments, while the relative abundances of Acidobacteria and Chloroflexi were higher in non-impacted sediment. High relative abundances of Bacteroidetes (class Bacteroidia) have been found to occur in flocculent matter and microbial mat samples within a hard-bottom salmonid farm in the Hermitage Bay area (Newfoundland, Canada), which at the time of sampling had been inactive for 3 months (Verhoeven et al., 2016). Flocculates sampled near salmonid farming cages exhibited higher abundances of Bacteroidetes than samples collected from less disturbed sites (Verhoeven et al., 2018). Since our community analyses focused on the surficial centimeter of the sediment and the samples typically consisted of fine mud, they likely included materials that had originally remained suspended in FIGURE 7 | Distance-based redundancy analyses (db-RDA) of sediment bacterial and eukaryotic communities, based on two types of models (see Table 4 for details). Data are shown for (A) bacterial communities (Model 1), (B) eukaryotic communities (Model 1), (C) bacterial communities (Model 2) and (D) eukaryotic communities (Model 2). Environmental variables were projected onto the ordination space using the R function envfit() fitted on linear scores of the ordination axes.
the water column (Milligan and Law, 2005). With the shifts in bacterial taxon composition detected in our study bearing many similarities to those reported in relation to active fish farming activities (despite our study site not being actively used for farming since 2008), the results point toward the possibility of a legacy impact of aquaculture-associated organic loading on the structure and composition of sediment bacterial communities.
In contrast with our results and those reported in Quero et al. (2020), Moncada et al. (2019) reported an increased abundance of Chloroflexi near fish farming cages. While we did not find evidence for a similar increase in the relative abundance of Chloroflexi, this could be due to site-specific factors. The decline we observed in the diversity of sediment metazoa near the farm was, however, in agreement with general evidence for a shift toward a microbially-dominated state within areas impacted by aquaculture (La Rosa et al., 2001;Mirto et al., 2012). Interestingly, members of the class Flavobacteriia were detected even in sites with comparatively high sediment depth-integrated O 2 consumption rates and low bottom water O 2 concentrations ( Table 3), despite this class often being regarded as aerobic (Kirchman, 2002;Kumagai et al., 2018). The physiological mechanisms enabling this class to persist in low-oxygen environments have been extensively characterized by Kumagai et al. (2018).
Our data showed distinct correlations between bacterial and eukaryotic community structure in relation to the sedimentological and water column variables that were examined (Tables 3, 6). The results further indicated that shifts in the structure of both community types were correlated to a shared subset of environmental variables, including sediment porewater

NH 4
+ concentrations and the sediment depth-integrated oxygen consumption rate. While these and other variables have been previously related to variation in either prokaryotic (bacterial and archaeal) or eukaryotic community structure within eutrophic sediments (Edlund et al., 2006;Chariton et al., 2015;Dowle et al., 2015;Fodelianakis et al., 2015), our cross-comparison of 16S and 18S rRNA metabarcoding data provides simultaneous information about both domains. The net intensity of organic matter remineralization, as recorded by the porewater NH 4 + inventory and rate of O 2 consumption, appears to be a key predictor of both bacterial and eukaryotic community structure in the studied locations. Intensity of remineralization is strongly coupled to microbial activity and hence the measured environmental variables.
Considered together, our data support the use of both 16S and 18S rRNA gene metabarcoding analyses as efficient tools to detect anthropogenic impacts such as eutrophication on marine ecosystems Stoeck et al., 2018). The overlap in bacterial and eukaryotic community responses to nutrient enrichment and oxygen depletion suggests that either method can provide useful insights into the overall impacts of eutrophication on ecosystem health. Where the variable of interest is known to be associated with shifts in the structure of both bacterial and eukaryotic communities (such as in the case of NH 4 + ), 16S rRNA metabarcoding can provide a highly sensitive measure of ecosystem-level responses to eutrophication. The bacterial community data also exhibited less within-station variation than the eukaryote community data (with significant variation in multivariate dispersions only observed for the eukaryote data set; Figure 6). While it remains possible that such differences are partially influenced by the use of small sediment volumes (per-replicate mass of 0.25 g; "Materials and Methods" section) (Nascimento et al., 2018), rarefaction analysis and Good's coverage estimates for OTUs were indicative of sufficient sampling depth in both the 16S and 18S rRNA gene data sets in our study (Supplementary Figure 1 and Results). Ultimately the choice of which metabarcoding approach to use will depend on the overall goals of the monitoring survey, as well as the available resources. Sediment grain size, for example, was found to have a stronger influence on eukaryotic community structure than on bacterial community structure (Table 6), which could be of relevance to interpreting the ecological impacts of activities that can disturb the surrounding sediment grain size composition, such as drilling and dredging (Smit et al., 2008) as well as offshore wind farms (Floeter et al., 2017). The sediment C org :N tot ratio was also more strongly associated with eukaryotic than bacterial community shifts, while variation in bacterial community structure was particularly tightly coupled to sediment porewater NH 4 + concentrations. The results of this study further suggest that the selection of environmental metadata variables to measure can have wideranging consequences for interpreting the results of marine impact assessments. For both the bacterial and eukaryotic data sets, db-RDAs using variables based on their low collinearity explained a greater proportion of the variance in community structure than the second set of db-RDAs employing alternative variables (Tables 4, 6). This implies that, although measurements performed using bottom water samples did correlate with the structure of sediment bacterial and eukaryotic communities, the sensitivity of environmental monitoring protocols could be improved by substituting these variables with appropriate sedimentological data (such as sediment porewater nutrient inventories and/or depth-integrated O 2 consumption rates). The practical feasibility of such modifications is expected to depend on the availability of field equipment, including sediment coring devices and oxygen microelectrodes, and the depth resolution of porewater sampling, which impacts on overall sample processing time. Analytical approaches for porewater nutrient determinations are comparable to those of bottom water samples and hence do not introduce new technical requirements in the laboratory.
The current study focused on a single field transect with samples collected during a single month (see "Materials and Methods" section). Therefore, to ascertain the generality of our findings and to identify sets of low-collinearity variables that could be efficiently used for monitoring purposes across a broad range of settings (including several types of habitat and different seasons), the conceptual approach presented herein could be extended to alternative locations and studies conducted over longer timescales. To further improve our ability to devise the best possible tools and strategies for benthic marine impact assessments, the types of comparisons performed in our study could also be applied to functional gene expression data. Moncada et al. (2019) showed a reduction in the abundance of dissimilatory sulfate reductase and nitrite reductase genes in mariculture-impacted sediments. While the authors were able to establish correlations between gene expression data and the abundances of selected microbial taxa, they noted that in certain instances this was impossible, likely because we only have a limited understanding of the ecological functions performed by many of the OTUs within sequencing data sets (Moncada et al., 2019). While rRNA gene metabarcoding analyses do not yield direct information on functional gene expression, the choice of which metabarcoding approach to use (and what environmental variables to measure) could nevertheless be greatly aided by the availability of such data.
In summary, our results show that rRNA gene metabarcoding analyses targeting bacterial or eukaryotic communities provide congruent information on the impacts of aquaculture-associated eutrophication on sediment communities, with shifts in the structure of both community types being correlated with a common set of sedimentological and water column variables. Ultimately the choice between 16S and 18S rRNA gene metabarcoding as a tool for monitoring the health of natural ecosystems should be guided by the aims of the study in question and the environmental variable(s) of interest. In this context, our findings clearly demonstrate how the usefulness of both metabarcoding methods is directly dependent on the types of environmental metadata that have been collected. To facilitate the adoption of eDNA metabarcoding as a tool for marine impact assessment, a major challenge will be to determine how diverse physical and chemical variables correlate with microbial versus macrobial community shifts not only locally, but regionally, seasonally and with reference to multiple habitats.

DATA AVAILABILITY STATEMENT
Raw 16S rRNA and 18S rRNA gene sequences supporting the results of this article are available in the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/docs/) under BioProject accession number PRJNA725053. R scripts used for data analysis, phyloseq objects including data processed in micca and environmental metadata (Materials and Methods) are available on GitHub (https://github.com/jessepharrison/ seili-metabarcoding). A Singularity container (https://sylabs.io/ singularity/; Kurtzer et al., 2017) definition file and package installation script for reproducing the R environment and package versions used for data analysis are provided in the repository.

AUTHOR CONTRIBUTIONS
JH performed DNA isolations and initial PCRs for sequencing analyses, bioinformatics and data analyses, and wrote the manuscript. KK was responsible for the study design. KK, IS, and P-MC carried out the fieldwork and sampling. KK and TJ performed sample preparations for geochemical analyses and performed the HS measurements. All authors contributed to data interpretation, reviewed the manuscript and approved of its publication.

FUNDING
This work was funded by the Academy of Finland (Grant Nos. 278827 and 304006 awarded to KK).

ACKNOWLEDGMENTS
We thank the R/V Aurelia crew, Markus Weckström, and Anni Jylhä-Vuorio for help with sample collection, and the staff of the Archipelago Research Institute (University of Turku, Finland) for access to sample processing facilities. We acknowledge Lars Paulin and the staff of the DNA Sequencing and Genomics Laboratory (Institute of Biotechnology, University of Helsinki, Finland) for their guidance and carrying out the Illumina MiSeq sequencing analyses. Helen de Waard from Utrecht Geolab (NL) is acknowledged for elemental analyses of the porewater samples and Tvärminne Zoological station for nutrient analyses of porewater samples.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021. 708716/full#supplementary-material Supplementary Figure 1 | Rarefaction curves for (A) bacterial and (B) eukaryotic Illumina MiSeq metabarcoding data. The data were generated following PCR amplification of bacterial 16S rRNA genes and eukaryotic 18S rRNA genes amplified from eDNA isolated from sediments from six sampling sites in the Archipelago Sea (Baltic Sea, Finland). Curves are shown for OTUs based on a similarity cut-off threshold of 97%, following removal of chimeric sequences.
Supplementary Figure 2 | Non-metric multidimensional scaling (nMDS) ordinations of sediment communities in field sites sampled along a eutrophication gradient in the Archipelago Sea (Baltic Sea, Finland). The ordinations were derived from Aitchison resemblance matrices calculated from Illumina MiSeq OTU abundance data. Data for bacterial communities are shown (A) prior to filtering out sequences assigned as NA, Chloroplast (class level) or Mitochondria (family level), (B) after filtering out NA sequences, and (C) after filtering out NA, chloroplast and mitochondrial sequences. Data for eukaryotic communities are shown (D) prior to filtering out NA sequences and (E) after filtering out NA sequences.