Environmental drivers of eukaryotic plankton and fish biodiversity in an Arctic fjord

Arctic fjords are experiencing rapid environmental shifts due to climate change, which may have significant impacts on marine biodiversity and ecosystem functioning. However, the impact of climate change on fjord biodiversity is difficult to quantify given the low accessibility and high cost to sample these areas. In this study, we sampled locations from inside to outside an Arctic fjord and used environmental DNA metabarcoding to assess how the biodiversity of fish and eukaryotic plankton communities relate with environmental conditions. We detected a total of 12 fish taxa and 872 MOTUs for eukaryotes and found marked gradients of temperature and salinity driven by the distance to the glacier at the terminal part of the fjord and depth. Eukaryotic richness was mainly associated with lower temperature and chlorophyll a. Moreover, co-inertia analyses showed a shared structure between eukaryotes assemblages and the environmental gradients. A partial association between eukaryote and fishes suggest a potential shared effect of environmental gradients among these taxa. Our findings provide a baseline for future studies to assess how these assemblages may be impacted by ongoing environmental changes and highlight how fjord ice loss and warming might shift environmental gradients and species distribution under climate change.


Introduction
In the coastal Arctic region, the connection between land and sea is characterized by the flow of meltwater and runoff into fjords, which are deep, glacially carved estuaries (Bianchi et al. 2020). Water flows either from surface streams or subglacial discharge at the calving front. Coastal glaciers in this region either terminate in fjords, known as marineterminating glaciers (or tidewater glaciers), or on land, leading to proglacial streams that drain into the fjord, referred to as land-terminating glaciers, driving distinct oceanographic circulations (Meire et al. 2017;Henson et al. 2023). Fjords display a seasonally dynamic, stratified circulation with a buoyant surface layer that flows offshore, compensated by the inshore advection of (subsurface) oceanic water (Haine et al. 2015). This cold water from melting ice causes marine water circulation funneling large fluxes of seasonal, buoyant freshwater offshore, influencing circulation dynamics (Haine et al. 2015), and shaping complex gradients of salinity and productivity, with important biomass generation at all trophic levels (Piwosz et al. 2009). As a result, fjords harbor diverse and abundant biodiversity, which serves important functions such as organic carbon burial (Smith et al. 2015). However, climate change is modifying the conditions of fjord systems with increased glacier melt and discharge as well as declining sea ice leading to a transformation of the marine ecosystem (Kwok et al. 2009;Smith et al. 2015;Meire et al. 2017;Holding et al. 2019;Hopwood et al. 2020;Hawkings et al. 2021). These ecosystem modifications due 1 3 to climate change may have profound effects on the fjord marine food web, and affect the availability of resources for fisheries and human livelihoods (Hansen et al. 2003). Therefore, understanding the dynamics of fjord ecosystems is critical for maintaining biodiversity, particularly as ongoing climate change is likely to further disrupt ecosystem functioning. (Babin 2020).
Fjord ecosystems in the Arctic are already impacted by ongoing climate change with shifting species composition (Willis et al. 2008;Li et al. 2009;Fujiwara et al. 2011;Kortsch et al. 2012), but those changes are difficult to monitor due to the lack of accessibility of these environments. Hence, traditional monitoring survey methods cannot keep pace with the rapid environmental changes occurring, and the rocky bottoms of fjords make it difficult to use approaches such as bottom trawling while bioacoustics shows limited taxonomic resolution (Stanton 2012). Alternative methods relying on images through underwater cameras or ROVs come with the advantage of being non-destructive, but the deep waters prevent most classic cameras to be used, whereas ROVs deployments are very expensive and can only cover a small spatial area per deployment. Environmental DNA (eDNA) metabarcoding has emerged as a powerful method to survey biodiversity, from micro to macro-organisms (Pawlowski et al. 2021;Takahashi et al. 2023). All living organisms shed DNA in their environment through the production of mucus, gametes, urine, excrement, and more (Deiner et al. 2016). Traces of this DNA can be found in the environment and can be detected from a small environmental sample (water, soil, etc.) generally up to a few hours up to a few days after it was shed in the water (Collins et al. 2018;Andruszkiewicz Allan et al. 2020;Holman et al. 2021). Note that environmental DNA samples are complex assemblages of DNA in many states (including entire organisms under a certain size class; see Pawlowski et al. 2020;Rodriguez-Ezpeleta et al. 2021). Because of the persistence of eDNA, just a few liters of water is often all that is needed to capture this eDNA. After filtration, the eDNA can be extracted, amplified, and sequenced using metabarcoding while being non-invasive for many macro-organisms. With bioinformatics tools, the DNA sequences can be sorted, cleaned, and compared to a reference genetic database to identify the species present in the environment (Boyer et al. 2016). There is now a considerable body of literature using eDNA metabarcoding from marine water samples in most ecosystems, from cold Arctic (Merten et al. 2023) and Antarctica (Fonseca et al. 2022) to temperate (Rozanski et al. 2022) and tropical environments (Polanco Fernández et al. 2020) targeting generally fish but also other groups such as eukaryotes, sharks, mammals. A few tree-of-life approaches have shown that we can recover a large taxonomic breadth of species from just a single water sample using a multi-primer approach (Stat et al. 2017). The eDNA field is now moving forward from proofs of concept toward applications for long-term monitoring of biodiversity (Ruppert et al. 2019).
Svalbard fjords are particularly affected by ongoing climatic changes and are transforming faster than most marine environments. Higher latitudes have warmed four times faster than the global average since 1979 (Rantanen et al. 2022), with local variation up to six times notably in Svalbard (Adakudlu et al. 2019). Svalbard fjords are classified as subpolar and freeze periodically during winter. Nowadays, sea ice decline has been widely documented around the archipelago even during winter (Onarheim et al. 2014) and numerous fjords almost do not freeze anymore due to the phenomenon of 'atlantification' (Kujawa et al. 2021). The western Svalbard shelf features numerous troughs, 200-400 m deep, allowing for water exchange between the open ocean and fjords, which rapidly change environmental conditions. Western Svalbard fjords have a complex circulation, with a dominated Arctic water current influence except during the summer months when it shifts to warmer and more dense water of the Atlantic water current. Global changes have modified water mass flows, and warm Atlantic water now often enters the fjord during the winter instead of cold Arctic waters (Cottier et al. 2007). This led to several ice-free winters and have influenced abiotic fluxes and the abundance and distribution of biological communities (Willis et al. 2008;Hop et al. 2019). To date, many studies come from the southern Kongsfjorden fjord due to the proximity of a research station (Wiencke and Hop 2016), but to be able to generalize the understanding of how complex circulation scenarios in fjords structure biological communities, we need to further investigate whether they hold in both similar as well as divergent habitats. Furthermore, many studies tend to focus on a single branch of the tree of life and rarely assess multiple trophic levels simultaneously. (Basedow et al. 2004;Bourgeois et al. 2016;Bhaskar et al. 2020). Due to these spatial and taxonomic biases, we have a relatively limited picture of how these changes have impacted their ecosystem from micro to macro-organisms. Metabarcoding of eDNA enables us to decipher patterns in community structure across divergent taxa of the tree of life from the same eDNA samples. To get a better understanding of the tangled links between the environmental conditions and biodiversity, it is necessary to gather basic knowledge on its biodiversity structure along environmental gradients within these fjords.
Here we investigated the relationship between environmental parameters and biological composition from plankton to high-trophic level fishes within an Arctic fjord system using eDNA metabarcoding. In this study, we sampled the Lilliehöök Arctic fjord in Svalbard (79 N) from the glacial front up to the fjord mouth from the surface to a maximum of 85 m depth. We sampled water for eDNA metabarcoding targeting planktonic and fish communities with two sets of universal primers. Alongside water samples, we measured abiotic parameters on the water column using a CTD and multiparameter probe. We asked the following research questions: (1) What are the main physico-chemical environmental gradients from the inside to the outside of the fjord and from the surface to the depth and does it relate to the distance from the glacier? (2) How are the environmental gradients associated with planktonic and fish assemblages within the Fjord? (3) How are fish and planktonic diversity structured along the fjord from the mouth of the fjord to the ice-terminating section?

Study area and sampling
We sampled the Lilliehöök fjord on the west coast of Spitsbergen (Svalbard, Norway) over 3 days from 3 to 5 of August 2021. Samples were taken from the glacier front up to the fjord mouth of the Krossfjorden system, around 30 km long, after the Lilliehöök fjord merged with the mouth of Möller fjord. The fjord's maximum depth has been recorded at 373 m (Svendsen et al. 2002) and has no sill at its entrance, thereby facilitating water exchange with the open ocean of the West Spitsbergen Current. We used a research vessel to sample 5 sites for a total of 15 samples, sampling 3 depths per site (3 m, chlorophyll a maximum and 85 m, unless sea floor was shallower). Shallow and intermediate samples between 3 and 12 m represent ~ 35 L of water filtered in situ using long tubing and a peristaltic pump, and all other deeper samples were taken from a total of 3 Niskin bottles (General Oceanics), representing 22 L of water sampled per sample. Water was filtered through a VigiDNA filtration capsule (SPYGEN) with a 0.20 µm pore size using an Athena peristaltic pump (Proactive Environmental Products, Bradenton, Florida) with a flow rate of ~ 1 L min −1 . Each sample was handled with single use tubing and gloves.

Analysis of physical parameters
The structure of the water column was characterized using a CTD (Idronaut Ocean 316) equipped with sensors for conductivity, temperature, PAR (Photosynthetic Active Radiation) (sensor xy), and chlorophyll a (chl a; sensor xy). The CTD was left 1 min at the surface to measure incident light out of the water before data acquisition. The conductivity profile for the section of the study site and the temperature salinity plot to identify changes in water masses were done using R by applying a multilevel B-splines method available in the MBA package from the methods of (Lee et al. 1997).

DNA extraction and amplification
We performed DNA extraction and amplification in distinct dedicated rooms set up with positive air pressure, UV treatment, and frequent air renewal. We processed the eDNA capsules at SPYGEN, following the protocol proposed by Polanco-Fernández et al. (2020). After DNA extraction, we diluted any samples identified as inhibited fivefold before amplification (Biggs et al. 2015)). We performed the DNA amplifications in a final volume of 25 μl, using 3 μl of DNA extract as the template. Teleo is located on the mitochondrial 12S region, designed to amplify both Actinopterygii and Chondrichthyes DNA (Polanco-Fernández et al. 2020). For teleo, we 5′-labeled the primers with an eight-nucleotide tag unique to each PCR replicate, allowing the assignment of each sequence to the corresponding sample. The tags for the forward and reverse primers were identical for each PCR replicate. We ran 12 PCR replicates per sample to account for rare species and stochasticity. We denatured the PCR mixture at 95 °C for 10 min, followed by 50 cycles of 30 s at 95 °C, 30 s at 55 °C, and 1 min at 72 °C, and we completed a final elongation step at 72 °C for 7 min. We pooled the purified PCR products into equal volumes to achieve a theoretical sequencing depth of 1,000,000 reads per sample. The 1389F/1510R primer pair, hereafter referred to as V9-18S, is located on the 18S gene and designed to amplify a large breadth of eukaryotic lineages. It is now a recognized standard to study eukaryotic biodiversity and is commonly used by long-term ecological series and worldwide plankton exploration (Amaral-Zettler et al. 2009;de Vargas et al. 2015). While this primer pair is designed to amplify eukaryotes, note that its large breadth of amplification can also amplify non-target groups such as some prokaryotes. We 5′-labeled the primers with an eightnucleotide tag unique to each sample, with identical tags in forward and reverse. We ran 12 PCR replicates per sample. We denatured the PCR mixture at 95 °C for 10 min, followed by 40 cycles of 30 s at 95 °C, 30 s at 57 °C, and 1 min at 72 °C, and we completed a final elongation step at 72 °C for 7 min. We pooled the purified PCR products into equal volumes to achieve a theoretical sequencing depth of 250,000 reads per sample. For V9-18S, we 5′-labeled the primers with an eight-nucleotide tag unique to each sample, allowing the assignment of each sequence to the corresponding sample. The tags for the forward and reverse primers were identical for each sample. After amplification, we quantified the samples using capillary electrophoresis (QIAxcel; Qiagen GmbH, Hilden, Germany), and we purified them using the MinElute PCR purification kit (Qiagen GmbH). Before sequencing, we quantified the purified DNA again using capillary electrophoresis. Two libraries (one for each marker) were prepared by DNAGensee (Le Bourget du Lac, France) using the TruSeq PCR-Free kit (Illumina) and sequenced them separately. We carried out paired-end sequencing using a MiSeq sequencer (2 × 125 bp, Illumina, San Diego, CA, USA) on the MiSeq Flow Cell Kits (v3; Illumina), following the manufacturer's instructions. We amplified (12 replicates) and sequenced two negative extraction controls and two negative PCR controls (ultrapure water, one for each marker) in parallel to the samples to monitor possible contamination.

Bioinformatics
The 18S-V9 assay targeting planktonic organisms was processed to generate a list of MOTUs to approximate their richness. We followed the processing pipeline of the TARA studies using the same marker (de Vargas et al. 2015) where they used the SWARM clustering algorithm (Mahé et al. 2015). The entire workflow consisted of aligning pairedend reads with vsearch (Rognes et al. 2016), then sample demultiplexing and primer trimming using cutadapt (no mismatch tolerated for demultiplexing and two mismatches tolerated for primer trimming; (Martin 2011). Then we clustered sequences using SWARM (d = 1) and detected chimeras using uchime (Edgar et al. 2011). We further curated reads by removing all occurrences below 3 reads (de Vargas et al. 2015). Taxonomic assignments were matched to the protist reference ribosomal databases (PR2; (Vaulot 2022), a curated genetic database for the 18S-v9 primers and using the RDP classifier algorithm from the dada2 pipeline with defaults settings (Callahan et al. 2016). We accounted for tag-jump by removing all intra-library occurrences of less than 0.001 ratio of the total sequence abundance. Taxonomic assignments were then checked to remove any non-target amplification such as Homo sapiens and all sequences occurings in the blanks were removed from the data set.
We also discarded all sequences not assigned to any domain (e.g., bacteria, archaea, or eukaryotes) and only retained sequences assigned to eukaryotes for analysis. The fish sequence data with the teleo primers was processed to generate a list of taxonomic entities (taxa) corresponding to distinct species. We used a pipeline using the obitools toolkit (Boyer et al. 2016) encapsulated in the snakemake pipeline manager (https:// gitlab. mbb. univ-montp2. fr/ edna/ snake make_ rapid run_ obito ols, version 1.2.0). The main steps include the merging of forward and reverse reads using the illuminapairedend function, demultiplexing using nsgfilter, and cleaning using obiclean with default values. Taxonomic assignment was done using ecotag on the EMBL snapshot of May 2021. We investigated all perfect matches corresponding to several species to attempt getting finer taxonomic assignments, if for example, 2 species of the same genus corresponded to a sequence detected by only was known to occur in the Arctic, we replaced the genus-level assignment to the species-level by the species occurring locally. We accounted for tag-jump by removing all sequence occurrences below 0.001 ratio of their total abundance in a given library (Schnell et al. 2015). Data was further cleaned by removing all occurrences below 10-reads for a given PCR replicate, and we cleaned taxonomic assignments for taxonomic redundancy to better approximate species identity and richness under an incomplete genetic reference database. So, for each filter, any detection at a higher taxonomic level that was not a 100% perfect match was removed if a detection was present at a lower level. For example, in a given sample, we detected the species Clupea harengus (100% match) and the genus Clupea (98% match) so we removed the Clupea occurrence to have a representative number of taxa, as the genus detection is likely an error derived from the correct species-level sequence. This approach allows to get a conservative estimate of the number of species, considering the low fish richness in these latitudes, both overall and per clade, together with the coverage of the reference database (48% at least, as shown in GAPeDNA, (Marques et al. 2021). In cases with lower reference database coverage and higher richness, there are good reasons to use clustering, for example to better approximate richness, but in this particular case using a non-clustering pipeline was more appropriate.

Correspondence between environmental conditions and species composition
We used several analytical methods to examine the relationship between environmental parameters and taxa/MOTUs composition in eDNA samples. First, environmental parameter values for each eDNA sample were ordinated with a principal component analysis (PCA) using the FactomineR R package. From the biological compositional data, the taxa and MOTUs composition were ordinated using principal coordinate analysis (PCoA) using the ape R package. We also ordinated the decomposition of the beta-diversity by its total, turnover and nestedness components using betapart (Baselga and Orme 2012). Then, a co-inertia analysis was conducted using the coinertia function of the ade4 R package to evaluate variance variation between (i) fish and environment, (ii) eukaryotes and environment, and (iii) fish and eukaryotes. To determine the significance of the RV coefficient obtained, a Monte-Carlo test was performed on the sum of eigenvalues with 9999 permutations using the RV.rtest function.

Taxa richness gradient from inside to outside the fjord
We analyzed the relationship between local taxa richness and environmental parameters for plankton and fish. We computed local taxa richness as the number of MOTUs for plankton and the number of taxa for fish. We assessed the effects of standardized environmental parameters on local taxa richness using general linear model (GLM) with a poisson distribution and log-link function. Explanatory variables were checked for collinearity by computing correlation and checking for variable inflation factor (VIF) from the cars R package. Pairs of variables with a VIF superior to five were investigated to only keep one to limit the effects of multicollinearity.

Environmental gradients from inside to outside
The spatial gradient of environment parameters revealed a significant contrast between the sites F1-F2 and the rest of the samples due to the glacier's influence (Fig. 1), resulting in lower temperatures in surface waters. Additionally, a station with distinct environmental conditions was observed in F4, where warmer temperatures of ~ 4 °C were recorded at depths ranging from ~ 60 to ~ 90 m when the surrounding waters were ~ 3 °C. Surface waters close to the glacier were found to be more saline compared to surface waters farther away, indicating a likely local upwelling. Light penetration in the water column decreases away from the glacier, from 653 μmol photons m −2 s −1 at F1 station to 75 μmol photons m −2 s −1 at F5 station at 3 m. Important chl a concentrations (> 5 µg L −1 ) were found only above 30 m, with a peak observed around 5 to 20 m, and lower chl a levels recorded near the glacier. The deep chlorophyll maximum depth increases away from the glacier, with a depth of 13 m at F1, 8 m at F3, and 6 m at F5. We show that salinity and temperature exhibit opposite trends but are highly correlated (Fig. 2). Samples are arranged according to depth and distance to the glacier, but not clearly clustering together. Note that F1 at 3 m depth and F2 at 13 m are poorly represented in the multidimensional space (cos 2 < 0.5). The combination of high salinity and low temperature was observed across all sites and depths, except for the surface. The deeper samples (85 m) primarily clustered around the high salinity and low temperature position of the space, while intermediate samples at shallow depths (F3-F4) showed more clustering in the high temperature and low salinity space.

Biodiversity recovered from eDNA metabarcoding
Following all bioinformatic cleaning, we retrieved a total of 1′810′073 reads for the teleo assay and 1′759′293 reads for the eukaryote assay. 12 cleaned fish taxa were detected with teleo (Online Resource 1) and 872 MOTUs for the eukaryote assay (Online Resource 1). Identified chimeras represented 72 sequences for eukaryotes, after all other cleaning steps but before discarding non-target sequences.
For fish, no species occurred in all 5 sites (Online Resource 2, Fig. S1) and typical Arctic species were detected such as Lycodes esmarkii, a bathydemersal species, Clupea harengus and Mallotus villosus, both benthopelagic species. Some taxa are detected mostly close to the glacier (F1, F2) like Stichaeidae, Lycodes esmarkii or Leptoclinus maculatus, other are detected in intermediate or outer fjord area like Mallotus villosus, Liparis sp. or Amblyraja sp. and a few are detected in all areas like Cottidae, Gadidae or Cyclopterus lumpus. Mean beta diversity is high, with 0.79 for total, 0.69 for turnover and 0.10 for nestedness components, so it is mostly driven by taxa turnover component, representing species replacement.
For eukaryotes, about half of all division occurred in every sample and depth (Fig. 3). Others such as Foraminifera were detected only in F1 and F2 at 35 and 85 m depth, respectively. In terms of read counts, Metazoa (encompassing mainly copepod zooplankton) is the most abundant group for all samples except in shallower waters of sites away from the glacier (F4 and F5; Online Resource 2, Fig. S2), but not in terms of MOTUs richness locally. The Ocrophyta division, comprising photosynthetic algae, is abundant only in shallow waters in all sites (max 6 m depth). Another abundant group is the Dinoflagellata, present everywhere and abundant in all sites and all depths except the most shallow and representing the highest MOTUs richness locally in the deeper samples. Mean beta diversity is high (0.60), with 0.48 for MOTU turnover and 0.11 for nestedness. For both assays, accumulation curves show that sampling was sufficient to recover most of the biodiversity locally (Online Resource 2, Fig. S3).

Correspondence between environmental conditions and species composition
Eukaryote composition is ordinated by both distance to glacier and depth, first by depth where all deeper samples cluster together and then by distance to glacier (Fig. 4). PCoA x-axis is driven by temperature (pearson's correlation, r = − 0.89), salinity (r = 0.97), Chl a (r = − 0.55), depth (r = 0.91), PAR (r = − 0.43), whereas PCoA y-axis is mostly driven by distance to glacier (r = − 0.65), and temperature (r = − 0.23).
The site F4 at 85-m depth is an outlier and investigation of beta-diversity indicated that it was due to low diversity and driven by nestedness rather than MOTU turnover, so a shared set of MOTUs rather than a composition exclusive to this sample. Fish composition presents little pattern and few   samples, mostly clustered between close to the glacier (F1, F2) and further away (F3, F4). F5 has little richness so is likely not fully representative of the local community. PCoA x-axis is mostly driven by temperature (pearson's correlation, r = 0.81), salinity (r = − 0.84), and PAR (r = − 0.37), whereas PCoA axis y is mostly driven by chl a (r = 0.77), distance to the glacier (r = 0.47), and temperature (r = 0.53).
The co-inertia analysis of variance between biological communities and environmental parameters revealed that the co-inertia between fish and eukaryotic communities seems important (RV = 0.67; p-value = 0.35) as well as the co-inertia between fish and environment (RV = 0.62; p-value = 0.23) but not significant. However, the co-inertia between planktonic communities and environmental parameters is high and significant (RV = 0.80; p-value = 0.0001) since it has a larger number of samples considered in the analysis.

Richness gradient from inside to outside the fjord and with depth
From the total number of eukaryotic MOTUs (872), the mean local MOTU richness was 302 MOTUs (± 86 sd).
Locally, eukaryote richness generally increased with depth, except for F4, which decreased at a depth of 85 m (Fig. 5, Online Resource 2, Fig. S4). Other than depth, there is no clear tendency for MOTU richness spatially across the fjord. Environmental predictors for the GLMs were selected to keep variables with a VIF < 5, thereby removing salinity and phycoerythrin for all models. Eukaryotic plankton richness seems to be negatively driven by light (r = − 0.13, p < 0.0001), temperature (r = − 0.082, p = 0.00046), and by chl a (r = − 0.079, p 0.00095), but not driven by depth (p = 0.66) (Online Resource 3). As we expect different plankton groups to be influenced by different mechanisms, we also modeled some groups recognized as important in Arctic plankton. We modeled the richness of Crustacea, comprising mostly zooplankton calanoids, Dinoflagellates and the phytoplanktonic group Ochrophyta with the same set of predictors. Crustacea richness was influenced negatively by temperature (r = − 0.25, p = 0.034) and light (r = − 0.33, p = 0.019), Dinoflagellates richness was only influenced by temperature (r = − 0.28, p < 0.0001) and light (r = − 0.11, p = 0.023) and Ocrophyta was only influenced by the depth (r = − 0.23, p = 0.0092). For the fish, we only analyzed the surface samples. Fish richness was generally low (total n = 13 taxa), with only two taxa found in F5, the site further away from the glacier and a mean richness of 4.6 (± 2.0 sd). At the junction between the two fjords (F3-F4) we detect a higher fish diversity (6 to 8 taxa compared to 2-4 elsewhere).

Discussion
The Arctic fjords in Svalbard provide unique habitats for marine biodiversity due to their complex oceanographic processes and the presence of glaciers fostering a unique productivity (Bianchi et al. 2020). In this study, we investigated the composition of marine biodiversity along the main environmental gradients in an Arctic fjord in Svalbard. Our results indicate that the spatial gradient of environmental parameters, specifically temperature, salinity, and oxygen, had a significant effect on the distribution of marine eukaryotes. In higher trophic levels, fishes also showed spatial organization with a peak of diversity in the intermediate position of the fjord. To our knowledge, quantitative analysis comparing fjords with land terminating and marine terminating glacier in similar regional settings was not done using eDNA together with key environmental parameters. Therefore, we place our findings in the context of an Arctic fjord with a marine terminating glacier. The tidewater glacier seemed to have an influence on the environmental conditions along the fjord, with lower temperatures and higher salinity values observed in surface waters close to the glacier. As the extent of glaciers is decreasing under climate change (Rignot and Thomas 2002;Meier et al. 2003), ice decline and glacier retreat on land are likely to drastically change the circulation and composition of fjords in the future.
The fjord system was structured by marked environmental gradients especially as regards to productivity and temperature, which have been observed in other Arctic fjord systems in summer with a tidewater glacier (Weydmann et al. 2013;Ershova et al. 2021;De Rovere et al. 2022). In particular, chl a concentrations were not restricted to a few meters but extended to ~ 25 m, with a peak observed around 5-20 m, a common feature of fjords with marine terminating glaciers in summer. Such local upwelling can be induced by subglacial freshwater discharge (Meire et al. 2017;Bianchi et al. 2020), but also due to freshwater input from the nearby small lake and local heating due to sunny days at time of sampling. These upwellings bring deep, nutrient-rich waters to the surface, thereby fuelling productivity (Meire et al. 2017;Kanna et al. 2022). As freshwater input observed from the nearby lake into the bay between the F2 and F3 sites seemed low with persistent good weather and in the absence of data on flow rate from the lake, we assumed that subglacial melt associated with the glacier was a major driver for the observed upwelling. Additionally, lower chl a concentrations are recorded near the glacier, which corresponds with documented influence of ice water melt structuring the fjord environmental gradients (Arendt et al. 2016). The studied fjord had relatively warm (2.5-6.3 °C) water flowing during the sampling summer month in August, indicating the expected seasonal shift from Arctic water influence in the winter months to Atlantic water influence during the summer (Cottier et al. 2007;David and Krishnan 2017). Water exchanges with the ocean are facilitated due to the absence of a sill at its entrance. The combination of high salinity and lower temperature (< 4 °C) was observed across all sites and depths, except for the surface. Arctic fjords waters are generally highly stratified during the summer months, with the few surface meters exhibiting higher temperature and lower salinity due to ice melt or river discharge (Sciascia et al. 2013). Sampling stations closer to the tidewater glacier indicate a local upwelling, as water masses are typical of deeper water masses with colder, more saline and less oxygenated water, showing less stratification (Kanna et al. 2022;Slater et al. 2022). Eukaryote compositional data revealed a spatial and depth structure in the eDNA composition. The communities found in deeper samples, where waters present higher salinity and lower temperature values, are distinct from communities in intermediate and shallow waters, characterized by higher temperature and lower salinity values. Our findings based locally measured parameters and eDNA community composition add to the body of literature showing that fjords with marine terminating glacier are structured by complex oceanographic processes, generate local upwelling and are associated with clear oceanographic structuring mechanisms for biological communities along the fjord (Svendsen et al. 2002). Marked oceanographic gradients should be associated with clear species turnover in marine organisms as observed for other local marine ecosystems (Aschan et al. 2013;Porter et al. 2017;Pecuchet et al. 2018). West et al. (2020) have shown that the composition signal of eDNA was largely spatially organized within an ecosystem, associated with environmental conditions (Jeunen et al. 2019;Cantera et al. 2022). Following those previous studies, we demonstrated that there were marked differences in eukaryotes and fish composition at the scale of the fjord with a mean turnover of 0.69 and 0.48 for fish and eukaryotes respectively, between neighboring sites distant by 7 to 11 km spatially and 5 to 80 m in depth for the eukaryotes. Our findings show that the eDNA signal is localized in association with the local environmental conditions, even in a dynamic and cold fjord system, as highlighted by recent research in the sub-Arctic (Guri et al. 2023). The co-inertia revealed an association between the eukaryotic planktonic communities and environmental parameters, and the GLM model confirmed that these communities shift mostly along the temperature gradient. These findings are consistent with previous studies that have demonstrated the importance of glacier-induced parameters gradients such as temperature and salinity (Fujiwara et al. 2011;Arimitsu et al. 2012). Other parameters unassessed in the present study such as nutrients like nitrates were commonly recognized as limiting in fjord systems and therefore strongly influence phytoplankton bloom and productivity (Hopwood et al. 2018;Halbach et al. 2019). While the co-inertia analysis between fish and eukaryote showed matching structure, the limited number of samples did not allow to evaluate the strength of the signal but literature suggest strong links between environmental conditions, productivity, and fish biodiversity in Arctic fjords (Meire et al. 2017). These findings highlight the importance of considering multiple factors, both biotic and abiotic, when investigating the drivers of community structure in aquatic systems (Tremblay et al. 2015).
The marked environmental gradient from the tidewater glacier to the mouth of the fjord is expected to be associated with gradient of diversity. We find that eukaryotes composition is structured by depth in the fjord, whereas the signal is less clear for MOTU richness perhaps due to some nonlinear effects alongside the depth gradient. Specifically, photosynthetic planktonic groups such as Ocrophyta were more abundant and speciose in the upper water layer, as expected considering their dependence on both light and nutrients, none of our environmental parameters explained their MOTU richness besides depth. We did not quantify nutrients, but we could expect a MOTU richness or abundance structuration by a combination of nitrates or phosphates with light in an Arctic fjord (Halbach et al. 2019;Randelhoff et al. 2020). The lack of relationship between light and Ocrophyta richness could be explained either from a conjoint requirement of sufficient light and nutrients, for which we are lacking values, or that locally biomass could be influenced by light alone but higher biomass do not necessarily correlate with higher richness, sometimes the opposite in case of a group-specific bloom. Arthropod zooplankton richness was mostly influenced by chl a, through a bottom-up mechanism as previous studies have shown (Kwasniewski et al. 2013). Dinoflagellate harbored more diversity and read count in deeper water layers, but this group of protozoans comprises 1 3 both auto and heterotrophic protists. They are known to be present during several seasons in Western Svalbard and their diversity and abundance seems to be driven by a mix of bottom-up, top-down and environmental variables such as temperature (Levinsen and Nielsen 2002;Sherr et al. 2003;Bhaskar et al. 2020). Hence, the link we observe between eukaryotic plankton and the environment can be indirect for some primary consumers and rather depend on phytoplankton. Plankton studies often quantify cell counts rather than species diversity, making a direct comparison with our work difficult as read counts out of metabarcoding cannot be used as a true proxy of biomass due to the method's bias (Elbrecht and Leese 2015; Lamb et al. 2019). Over all eukaryotic MOTUs, pairwise site dissimilarity was mostly driven by nestedness, not MOTU turnover, indicating environmental filter on communities rather than the co-existence of completely distinct assemblages in distinct environmental conditions. Our findings are consistent with previous research on the relationship between plankton and environmental factors in polar regions, with a primary producer community structured by glacier-influenced hydrographic parameters like temperature and higher trophic levels dependent on productivity more than environmental parameters (Kwasniewski et al. 2013;Kubiszyn et al. 2014;Halbach et al. 2019). Overall, our study contributes to the growing body of literature on biodiversity in polar regions and highlights the importance of understanding the underlying drivers of biodiversity in these fragile ecosystems.
A major limitation of the present study is that the sampling method had a limited spatial coverage of the fjord system. Although we selected five sampling sites within a 30 km long fjord system, the spatial resolution of the study is still limited. Additionally, the limited number of samples (15) collected from the five sampling sites may not accurately capture the biodiversity of the area, as there could be species that were not captured due to stochasticity but this low number of samples is in line with other small scale studies (e.g., Hegseth and Tverberg 2013;Zhang et al. 2019). This could affect the accuracy and representativeness of the results obtained in the study. Further, using eDNA metabarcoding allows a complete screening of the tree of life locally, but quantitative measures of abundance remain challenging to obtain, even though recent work have suggested new metrics to use read counts as a proxy to local abundance with less bias (Gong and Marchetti 2019;Martin et al. 2022;Guri et al. 2023). Given the promising results shown here, future studies should consider more fjords and more extended sampling locations within each fjord to get more generality and precision on the findings presented in the present study.
In conclusion, our study demonstrates the ability of eDNA metabarcoding to provide a comprehensive assessment of the biodiversity of fish and plankton communities in a Svalbard fjord. We demonstrate the marked association between environmental gradients and primary productivity with plankton and fish assemblages at a limited spatial scale, which is contributing to the high overall physical and biological diversity found in fjord systems. As the Arctic continues to warm, some tidewater glaciers, notably in Svalbard (Torsvik et al. 2019) are expected to retreat up to the land. Shifting to a land terminating glacier will thereby profoundly change the water circulation and biogeochemistry of these fjords with a drastic projected impact on biological communities due to a sharp decline in productivity (Meire et al. 2017). Associated changes in temperature and nutrient availability may impact the richness and composition of plankton communities, which in turn may have cascading effects on the fish community and the wider ecosystem up to bird populations (Grémillet et al. 2015;Cauvy-Fraunié and Dangles 2019). Our study provides a baseline for future studies to assess how these communities may be impacted by ongoing environmental changes. Overall, this work highlights the importance of understanding the environmental gradients that influence marine biodiversity in Arctic fjords given that those are playing a key role in structuring biological assemblages.
Acknowledgements The authors are grateful to the crew onboard the ship for help during the sampling.
Author contribution VM, LP, CA and KD conceived and designed research. EM and CH conducted the field work. AV did the molecular analysis. VM did the bioinformatics analysis. VM, CA and LP did the analysis. VM led the writing with inputs from all authors.
Funding Open Access funding provided by Lib4RI -Library for the Research Institutes within the ETH Domain: Eawag, Empa, PSI & WSL.

Declarations
Competing interests The authors do not declare any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.