Environmental DNA metabarcoding reveals seasonal and spatial variation in the vertebrate fauna of Ilulissat Icefjord, Greenland

: Ilulissat Icefjord in Greenland is experiencing the effects of climate change, with the Sermeq Kujalleq glacier being one of the fastest-moving and most productive ice streams in Greenland. This is likely affecting the distribution of species in the fjord, including those important to local fisheries. Due to heavy ice conditions, few studies on environmental and ecological conditions exist from the fjord. However, new techniques such as environmental DNA (eDNA) meta barcoding now allow deeper insight into the fjord system. Here, we combine local ecological knowledge with data on hydrographic conditions, stable isotopes ( δ 18 O), and eDNA metabarcoding to investigate the spatial and seasonal distribution of marine fish and mammals inside Ilulissat Icefjord. Our eDNA results support local observations that Arctic char migrate to the southern fjord during summer, harp seals forage in large herds in the fjord system, polar cod is the dominant prey fish in the area, and Greenland shark likely does not reside in the fjord system. Lower predation pressure in the Icefjord, due to the absence of Greenland shark and polar bears as well as limited fishing/hunting, is presumably one of the reasons why ringed seals and Greenland hal-ibut are larger in the Icefjord. Furthermore, our results indicate that in summer, the southern branch of the fjord system has a more diverse community of vertebrates and different water masses than the northern branch and main fjord, indicating a time lag between inflows to the different branches of the fjord system. Our approach highlights the value of combining local ecological knowledge with scientific research and represents a potential starting point for monitoring biological responses in Ilulissat Icefjord associated with climate-induced changes.


INTRODUCTION
Fjords with marine-terminating glaciers sustain high marine biomass productivity (Meire et al. 2017, Hopwood et al. 2020, as these glaciers induce mixing of the water masses through sub-glacial freshwater discharge, calving icebergs, and melting events. These events create an estuary-like circulation of the water masses, which brings nutrients up to the photic zone of the water column (Mortensen et al. , 2020. The meltwater from the glaciers contains essential minerals (Hawkings et al. 2016, Meire et al. 2016a, which enhances primary production (Rysgaard et al. 2012, Juul-Pedersen et al. 2015, Meire et al. 2016b) and leads to hotspots for marine biodiversity and high species abundance (Lydersen et al. 2014). Ilulissat Icefjord (previously known as Jakobshavn Icefjord, henceforth referred to as the Icefjord) has one of the most active marine-terminating glaciers in the Northern Hemisphere, namely the Sermeq Kujalleq glacier in the main fjord branch (Tang et al. 2004). This particular glacier has become a symbol of the global warming crisis, and the area was designated as a UNESCO World Heritage site in 2004. Despite this attention, the Icefjord is understudied; the most recent publication related to our research was conducted in Disco Bay, adjacent to the fjord. The water masses in Disco Bay have been shown to have high primary production levels (Jensen et al. 1999), high concentrations of the Nordic shrimp Pandalus borealis (Folmer & Pennington 2000), and high biomass of Greenland halibut Reinhardtius hippoglossoides (Wieland et al. 2007). Onethird of all Greenland halibut caught in the Disco Bay area are caught within and just outside the Icefjord (Nygaard 2019). Both Greenland halibut and ringed seals generally have a higher body mass inside the Icefjord compared to those caught in Disco Bay (Schiøtt et al. 2021), suggesting especially high productivity and growth in this area. Despite the importance of the Icefjord for fisheries and wildlife, and the rapid changes currently induced by climate change, few studies have investigated environmental and ecological conditions inside the fjord.
The Icefjord is closed off by large icebergs and sea ice throughout most of the year. These conditions have made surveys from research vessels (e.g. bottom trawl surveys) inside the fjord system virtually impossible. Hence, environmental DNA (eDNA) analysis of water samples collected inside the Icefjord represents an attractive alternative method. The eDNA approach can be used to non-invasively monitor biodiversity (Taberlet et al. 2018, Sigsgaard et al. 2020, Afzali et al. 2021, including at fine temporal ) and spatial scales (Port et al. 2016, Berger et al. 2020. Organisms leave DNA traces of their presence through shed skin cells, mucus, or feces that can be detected in their environment (Dejean et a. 2011, Taberlet et al. 2012. eDNA is degraded in the aquatic environment within days to weeks (Dejean et al. 2011, Thomsen et al. 2012, Brys et al. 2021). This fast degradation means that eDNA from water samples can be assumed to represent a relatively narrow time period of species presence. Indeed, previous studies have shown that eDNA from seawater samples reflects temporal changes in community compo-sition across seasons (e.g. Sigsgaard et al. 2017) and even between night and day ). In addition to being degraded, eDNA expelled into the sea is continuously being diluted, and eDNA compositions have been found to reflect local habitat types at scales of down to ~60 m (Port et al. 2016).
The main goal of this study was to investigate which fish and mammal species reside in the Icefjord and whether differences in the vertebrate species composition exist between different sampling locations and seasons. An interdisciplinary approach was chosen, which combined (1) eDNA analysis, (2) hydrographic data (CTD), (3) stable isotope analysis of water samples (δ 18 O), and (4) a parallel interview survey that targeted locals who use the Icefjord as a fishing and/or hunting area. The interview survey gave some insights into seasonal variations and changes that the locals have observed, likely as the result of climate change, which gave us expectations for which species we would find in the metabarcoding results and what seasonal changes to expect. This combined data enabled us to describe the seasonal and spatial variation of the marine ecosystem in the Icefjord in greater detail than with eDNA alone.

Study site
The Sermeq Kujalleq glacier in the main fjord of the Icefjord produces many large icebergs, which usually makes it nearly impossible to safely navigate these waters. The main fjord of the Icefjord, Kangia, has 2 main branches: Sikuiuitsoq to the north and Tasiussaq to the south (Fig. 1), and it has 4 marineterminating glaciers. The main fjord is the deepest part of the Icefjord with a uniform depth of 750−800 m, while the northern branch is 500 m at the deepest, and the southern branch is around 200 m at the deepest (Stevens et al. 2016). There is a sill at the entrance to the main fjord with a variable depth of 150−245 m, deeper at the northern part of the entrance (Gladish et al. 2015). At the entrance to the southern branch, another sill is present, which is situated at a depth of 70 m (Stevens et al. 2016). In winter, this shallow sill near the entrance to the southern branch creates strong tidal currents that prevent sea ice from forming, resulting in an open water polynya. The polynya roughly spans 6−25 km 2 in area (estimated from Sentinel 1 satellite images over the last 6 yr, as Sentinel 1 began imaging the area in July 2015, and calculated with the open-source geographical information system QGIS).

Local informants
In parallel with the eDNA sampling, an interview survey was conducted with local fishers and hunters in which 33 informants participated. Selected results from this survey -specifically knowledge about the fish and marine mammal community -are included here for comparison with the results from eDNA metabarcoding. Detailed descriptions about the interview setup, their observations of changes in biological and environmental conditions, and the adaptation of the locals to these changes can be found in Schiøtt et al. (2021).

CTD measurements, stable isotope analysis, and eDNA sampling
Water samples for stable isotopes, eDNA, and CTD measurements were obtained from the same location and time. A 2 l standard water sampler (Ruttner 2 l; Hydrobios) was used to collect water samples. eDNA samples were collected at 5 and 100 m depth, while samples for stable isotopes were collected at 5, 25, 50, 75, and 100 m depth. CTD measurements were done using a CastAway-CTD (SonTek US patent #8272262) that made continuous measurements from the surface down to 100 m. Detailed description of the sampling procedure can be found in Text S1−S3 in the Supplement at www.int-res.com/articles/suppl/m706p091_ supp.pdf. . White triangles in the main fjord indicate the presence of large icebergs. Location 1 is described in the main text as being the 'northern branch' of the fjord system, while location 2 is in the 'southern branch', location 3 as the 'main/inner' parts of the fjord, and location 4 as the 'mouth' of the fjord system. The location of the larger map is shown as a red square on the map of Greenland. Map made with the open geographical tool QGIS

Environmental DNA extraction and metabarcoding
DNA extractions were performed using the Qiagen DNeasy Blood & Tissue kit (spin column protocol with alterations according to Sigsgaard et al. 2020). DNA extractions were done in a laboratory dedicated to eDNA work in which decontamination routines were followed, including UV lights, wiping all surfaces with DNA AWAY (Thermo Scientific TM ) and subsequently with 70% ethanol before working with samples, and using separate buildings for pre-and post-PCR work. DNA extraction blanks and PCR blanks were included in the laboratory setup to control for contamination during handling.
Samples were PCR-amplified using the forward primers Elas02_F (5'-GTT GGT HAA TCT CGT GCC AGC-3') and Tele02_F (5'-AAA CTC GTG CCA GCC ACC-3') and the reverse primers Elas02_R (5'-CAT AGT AGG GTA TCT AAT CCT AGT TTG-3') and Tele02_R (5'-GGG TAT CTA ATC CCA GTT TG-3') (Taberlet et al. 2018), which target the mitochon drial 12S rRNA gene of elasmobranchs and bony fishes, respectively. Primers were twintagged using 6 nucleotide (nt) tags preceded by NN or NNN (De Barba et al. 2014). Each PCR run was performed using unique tags for each sample and each control. PCR reactions were performed using 10 μl HotStarTag MasterMix (Qiagen), 10 μl double-distilled H 2 O, 1 μl BSA (Bionordica), 1 μl primer mix, and 3 μl DNA extract. For the PCR, the following thermocycler conditions were used: 95°C for 15 min, 45 cycles of 94°C for 30 s, 50°C for 30 s, 72°C for 1 min, and a final elongation at 2°C for 5 min. Fragment sizes were verified with 2% agarose gel stained with Gelred TM . Four PCR replicates were performed for each sample using the same tag in each replicate. Four library pools were made, each containing one replicate of each sample and each blank control (4 μl of each sample/ control). The pools were purified using Qiagen's MinElute PCR purification kit, following the manufacturer's protocol. DNA con centration was measured on a Qubit fluorometer (Invitrogen) to make sure there was enough DNA in each library. Libraries were built using the Illumina Truseq DNA PCR-free Low Throughput kit. The libraries were sent to the Novogene laboratory in Cambridge, UK, where they were sequenced on an Illumina Nova Seq 6000 sequencer using 150 bp paired-end se quencing chemistry and requesting 5 GB of output per library.

Sequencing data analysis
We used the highly parallelized pipeline Meta -BarFlow  to efficiently process the eDNA sequencing data. The most recent version of this pipeline can be found at https://github. com/evaegelyng/MetaBarFlow; specific scripts used in the present study are available upon request. MetaBarFlow is built mainly on scripts from Frøslev et al. (2017) and the R package DADA2 (Callahan et al. 2016). Specifically, samples were demultiplexed using 'Cutadapt' (Martin 2011), trimmed with 'sickle' (Joshi & Fass 2011) and then processed with DADA2 to remove errors produced during PCR amplification and sequencing. The final sequences were searched against a locally downloaded version (released 6 November 2021) of the NCBI GenBank nt database (NCBI Resource Coordinators 2016) with the blastn algorithm (Altschul et al. 1990), specifying a maximum of 500 sequence hits, a minimum of 90% query coverage per high scoring segment pair, and a minimum of 80% sequence similarity. The sequences were then taxonomically classified using the R package 'taxizedb' (Chamberlain & Arendsee 2021), only including hits with 100% query coverage and similarity > 98%. Sequences were also searched against a local database, which included early access to reference sequences produced by Jensen et al. (2023) and sequences produced by Jensen et al. (2022). Taxonomic assignment of the sequences was then updated when pertinent. In cases of ambiguous best hits, a last common ancestor taxonomy was ascribed (e.g. Gadus morhua and G. macrocephalus were collectively listed as Gadus spp.). We then filtered out taxa that only occurred in 1 of 4 PCR replicates of a sample and taxa that were detected in higher relative read counts in a control sample than in any seawater sample. We also discarded 7 samples that failed to generate sufficient sequence data, but at least 4 samples were retained from each of the sampling events. Finally, we disregarded hits to terrestrial animals as well as birds, as these were not the target in our study.
Species rarefaction curves were performed on individual PCR reactions per sample to assess the sufficiency of sequencing depth, and species accumulation curves per sample were performed to assess the sufficiency of PCR replication. Individually se quenced PCR replicates were rarefied to the median number of reads across all PCR replicates, after which the 4 PCR replicates were aggregated into samples by summing the reads. Samples were then rarefied to the minimum number of reads found in any sample to limit the effect of differential sequencing depth, using the R package 'ROBITools' (v.0.1; LECA 2012). A heatmap was created using R packages 'pheatmap ' (v.1.0.12;Kolde 2018) and 'RColorBrewer' (v.1.1-2;Neuwirth 2014), and stacked barplots were plotted using the R package 'ggplot2' (v.3.2.1; Wickham 2016).

Partial distance-based redundancy analysis
Partial distance-based redundancy analysis (partial dbRDA) was carried out to evaluate the significance of environmental effects on taxa distribution (capscale function). Jaccard distances were calculated with the presence−absence community data, and the effects of water temperature, oxygen isotopes, salinity, depth, time, and location of sampling on taxa composition were tested. Prior to analysis, all continuous environmental variables (i.e. oxygen isotopes, temperature, and salinity) were standardized, whereas depth, time, and location of sampling were used as categorical variables. A single fish species, Boreogadus saida, was present at all sites and thus was removed from the analysis. Principal coordinates of neighbor matrices (PCNM function) were calculated based on the geographical coordinates of sites in order to partial-out inherent spatial effects. A permutational ANOVA-like test was performed to assess the significance of the partial dbRDA model, its axes, and environmental predictors after 999 permutations. All analyses were run with the R package 'vegan' v.2.5-7 (Oksanen et al. 2022).
Greenland halibut is the main fish targeted in the fishery in the Icefjord, while ringed seals and harp seals are the main mammal species caught in the Ice-fjord. Polar cod is targeted with nets and used as bait in the Greenland halibut fishery, while the remaining fish taxa are only caught as bycatch. It should be noted that Arctic char is only caught in the southern branch by fishers from Ilimanaq during the summer. Additional information regarding certain taxa for which the informants had detailed knowledge is found in Table S1.

CTD measurements
Lower temperatures (down to around −1.5°C) were generally observed in the surface layers than at greater depths (where temperatures were up to +1°C) in March/April 2019, December 2019, and February 2020 ( Fig. 2A). In July 2019, on the other hand, higher temperatures (up to + 8°C) were generally observed at the surface than at 100 m depth (up to +1°C). Similarly, summer measurements from Disco Bay showed higher temperatures in the upper 20 m in July 2019 ( Fig. 2A). At 100 m depth, the temperature remained close to +1°C at all stations and seasons, except in the main fjord during December 2019, where warmer water of up to + 2°C was ob served. Temperatures down to −1°C were observed around 40 m depth in the northern branch in July 2019, which was colder than the temperature at the other stations and seasons ( Fig. 2A). Salinity was relatively low in the surface layers at all stations and seasons and increased with depth ( Fig. 2B). In contrast, salinity was more or less stable with depth in Disco Bay during June and July 2019 (Fig. 2B). In addition, high salinity values (up to ~35 PSU) were also observed below 75 m depth in the southern branch in July 2019, higher than the values observed at other depths and stations (Fig. 2B). Temperature and salinity below ca. 80 m depth in the northern branch of the Icefjord resembled values observed in Disco Bay. The temperature−salinity (TS) properties of the water types ( Fig. 2C) showed that both the northern and southern branches had Coastal Water (CW) except in July 2019, when the southern branch had more saline water compared to the northern branch and Disco Bay, with salinity values close to 35 PSU. This water likely originated from Sub Polar Mode Water (SPMW) that had cooled, as SPMW has similar salinity values ).

Stable isotopes
At all locations, δ 18 O values decreased with decreasing salinity (

eDNA metabarcoding
A total of 520 M raw reads corresponding to 260 M read pairs were generated. Across the 12 libraries, we obtained a varying sequencing depth with 27−98 M reads library −1 (average of 43 M reads library −1 ). After filtration for contaminants and removal of terrestrial species and birds, we found a total of 38 taxa. These included 20 taxa of fish (conservative count excluding e.g. unspecific hits to family) and 7 species of marine mammals, all of which are known to occur in Greenland. Initial quantitative analyses revealed that the Elas02 and Tele02 primers yielded almost identical relative proportions of dominant taxa, so Fig. 4 only shows the results for the Tele02 primer. However, a few rare taxa were only found with one of the 2 primer sets, so presence−absence analyses were based on both Tele02 and Elas02 (Table 1, see Figs. 5 & 6). Note that quite a few of the detected taxa were also detected in our controls, but these were detected in lower read counts compared to seawater samples and were thus kept in the data (Tables S2−S7). Replicate samples resembled each other reasonably, with few exceptions as seen in Fig. S1.

Diversity of species and response to environmental variables
The lowest number of taxa (10) was observed in the inner part of the main fjord, sampled in December 2019, while the highest number of taxa (27) was observed in the southern branch in the samples from July 2019 (Table 1). The following taxa were observed in all or most of the samples: capelin Mallotus villosus, polar cod B. saida, cod Gadus spp., sculpins Myoxocephalus spp., Greenland halibut R. hippoglossoides, harp seal P. groenlandicus, and ringed seal P. hispida (Table 1). The remaining taxa were observed in only a few samples, in a specific field campaign, or at specific sampling depths. A full list of taxa can be seen in Table 1. Fig. 5 illustrates that a higher number of taxa were observed in the samples collected in the southern branch and that, while the first 15 taxa listed in the figure were observed in most samples, the remaining 23 were observed in only a few samples. Significant effects of location, time of sampling, depth, salinity, and temperature on eDNA composition were detected (p < 0.005), whereas oxygen isotopes was not significant. Partial dbRDA revealed that 30.6% of total eDNA variance was explained by these variables (i.e. constrained variables), while 15.2% was controlled by spatial structure and 54.2% was residual. Overall, communities were strongly affected by the time and location of sampling. Northern-and southern-fjord-associated communities showed higher similarity in February compared to July, during which higher water temperatures were recorded in the southern fjord.  Table 1). Gadus sequences could not be determined to species level, but higher relative read count contribution was particularly noticeable in the 5 m samples from the southern branch from April 2019 (Fig. 4). Some of the detected marine taxa showed seasonal differences in relative abundance. The most noticeable seasonal changes in respect of 3 taxa were as follows: (1) Arctic char S. alpinus had higher relative read counts in the samples from the southern branch from July 2019 at 5 m depth (Fig. 4); (2) harp seals P. groenlandicus had increased relative read counts in the samples from the northern and southern branch sampled in February 2020 (Fig. 4); (3) higher relative read counts were observed for the Greenland halibut R. hippoglossoides in the southern branch during July 2019 (Fig. 4). Whales generally showed low relative read count contributions (Fig. 4), and the common minke whale Balaenoptera acutorostrata, fin whale B. physalus, and narwhal Monodon monoceros were only observed in the samples from July 2019 (Table 1). However, the beluga Delphinapterus leucas showed high relative read counts in the samples collected from the main fjord (Fig. 4). Ringed seal P. hispida read counts were relatively lower than e.g. the harp seal, but somewhat higher relative read counts for ringed seals were observed in the southern branch compared to the northern branch and main fjord (Fig. 4), and this species was observed in all samples (Table 1)

DISCUSSION
Until now, ice conditions have prevented scientific field surveys of the marine ecosystem of the Icefjord, although locals from the interview survey have described changes related to season and climate change (Schiøtt et al. 2021). For example, locals report that (1) Arctic char migrate to the southern fjord during the summer; (2) large herds of harp seals come to feed on capelin in the main fjord; (3) Greenland halibut are larger in the main fjord of the Icefjord compared to Disco Bay and the northern and southern branch of the Icefjord; and (4) Greenland halibut are larger the further one moves into the main fjord, where the fish community is also less diverse (Table S1). Our study allows these reported differences to be scientifically investigated for the first time. Through eDNA metabarcoding, we were able to support the local knowledge of seasonal variations of Arctic char and harp seals and the lower diversity of taxa at the inner station in the main fjord. In addition, the CTD measurements and stable isotope analysis revealed different physical conditions between the branches of the fjord system, which likely contributed to the observed differences in fish communities.
Relative eDNA read count contributions of individual fish taxa have been found to correlate with relative abundance and biomass, both in controlled settings (Di Muri et al. 2020, Rourke et al. 2021) and in the field (Thomsen et al. 2016). We have therefore included relative read counts in our ecological analyses, but it is important to realize that several factors, like PCR bias, can influence relative read counts and thereby compromise inferences on relative biomass/abundance (Fonseca 2018, Yates et al. 2019. We also found that many of the same taxa observed in our field samples were also observed in the field controls as well as extraction controls. These could have been caused by a mix-up of primer tags, as some of the field sample replicates revealed relatively low eDNA concentrations (which, due to a mix-up of tags could have been field or extraction controls and not field samples), which lead to some of the samples only having 4 replicates instead of 5 (as seen in Fig. S1). However, as replicates corresponding to the same sample re semble each other, it can be assumed that a potential mixup of tags only happened to a few of the samples. Read counts in field and extraction controls were still lower than what was observed in the field samples, so this did not have a great impact on our overall eDNA results.

Distribution of species; spatial and temporal differences
Based on the results of the metabarcoding analysis in the present study, it is worth noticing that there are slightly different communities in each branch of the fjord system (Table 1, Figs. 4 & 5). The communities differed both in the presence−absence of taxa, but also in the relative read count abundance of the present taxa.
The most noticeable difference between the samples was in those collected in July 2019 in the southern branch and the northern branch as seen in Figs. 5 & 6, and also observed in the partial dbRDA that showed that the higher temperatures observed in the southern branch during July likely had an effect on the observed taxa (Figs. 6 & 7). Higher relative read counts are, for example, seen for Arctic char in the samples from the southern branch (Fig. 4), which can be explained by anadromous Arctic char spawning in ocean-connected lakes in the southern branch during June−July. Higher relative read counts in eDNA due to spawning events have previously been reported for marine fish (Tsuji & Shibata 2021), but in this study, the eDNA sampling was done at sea and not at the freshwater spawning site. However, the relatively higher read counts for Arctic char eDNA could also indicate a combination of spawning and migration. The Arctic char in Greenland and their complex variation in landlocked and anadromous lifecycles is poorly studied; this knowledge gap could potentially be addressed with eDNA in the future. Higher relative read counts for Greenland halibut eDNA are also noticeable in the July 2019 samples from the southern branch (Fig. 4), but as they are assumed to spawn during the period November−January (Albert et al. 2001, Greenland Institute of Natural Resources unpubl. data) the higher relative read counts were most likely not caused by spawning events, but rather by migration of Greenland halibut into the southern branch. Indeed, the entire stock of Greenland halibut in the Icefjord and other coastal Greenlandic areas results from the immigration of larvae and juveniles from offshore spawning sites in the Davis Strait (Vihtakari et al. 2022). In many Greenland fjords, the large adult fishes are thought to be non-spawning expatriates, 'trapped' in the fjords due to the shallow fjord entrance sills (Boje 1994). But contrary to this supposition, Boje et al. (2014) showed with acoustic tags that the Greenland halibut occupied warmer, shallow water in Disco Bay in summer but migrated into the Icefjord, to cooler and deeper water during winter and returned the fol-lowing spring. Locals have reported that Greenland halibut caught in the southern branch are noticeably smaller than the ones caught in the main fjord and that the fish community in the main fjord is less diverse, being dominated by larger Greenland halibut that prey on smaller fish including juveniles of 102  Fig. 7. Ordination of eDNA taxon scores according to partial distancebased redundancy analysis. Taxa are colored in black. Red arrows represent the estimated effect of salinity and temperature. First (RDA1) and second (RDA2) axes explained 45.9 and 19.7% of community data, respectively their own species. This led us to speculate that the southern branch might be a nursing ground for smaller Greenland halibut, taking refuge from predation. Greenland halibut generally show a very strong correlation between body size and depth, with larger fish inhabiting deeper waters (Bowering & Nedreaas 2000), so it would be expected that larger specimens are found mainly in the deeper main fjord and smaller ones are found in the shallower northern and southern branch. These higher relative read counts in Greenland halibut eDNA during July 2019 could suggest that juvenile Greenland halibut migrate into the southern branch during summer. An increase in the relative read counts for cod (Gadus spp.) was observed in the southern branch sampled in April 2019 (Fig. 4), which could not be identified to species level due to overlapping barcodes (equal score for G. morhua and G. macrocephalus). It is plausible that these higher relative read counts for cod are caused by spawning events by the Greenland cod, as they spawn around March− April (Mikhail & Welch 1989, Morin et al. 1991, which suggests that the southern branch is a nursing ground for Greenland cod as well. The temperature measured in the southern branch in April 2019 seems too low to support spawning events for the Atlantic cod ( Fig. 2A), as they require higher temperatures for the successful growth of their larvae (Laurence 1978, Otterlei et al. 1999. Altogether, this could suggest that the southern branch represents a nursing ground for Greenland cod in addition to Greenland halibut. Four informants described the fish community as being less diverse further into the main fjord, becoming dominated by large Greenland halibut that prey on smaller fish (Table S1). This information was supported by the eDNA results, as fewer taxa were represented in samples collected at the inner part of the main fjord compared to the other locations (Table 1, Fig. 5). The main fjord was not accessible during the other field campaigns, as it was closed off by massive icebergs and sea ice, so, unfortunately, this can not currently be confirmed. However, based on local knowledge, it can be assumed that the observed less diverse fish community in the inner part of the main fjord is the general pattern. The fishers also take advantage of the ice clearing events (when the main fjord clears of ice) and sail in as far as they can into the main fjord to catch larger Greenland halibut and thus achieve a higher fishing efficiency.
Interestingly, the locals have described that the deep-water fish roughhead grenadier has decreased drastically in abundance in recent years (Table S1), which could explain why this species was not ob -served in the metabarcoding results (Table 1, Fig. 5). However, since our water samples were collected at 5 and 100 m depth, taxa found at greater depths might not have been fully represented. However, as subglacial discharge from marine-terminating glaciers will entrain deeper water to the surface due to mixing , it can be assumed that fish that reside at greater depths would be represented in the eDNA metabarcoding results to some degree. There is some support for this concept; e.g. the presence of Amblyraja hyperborea (Table 1, Fig. 5).
The high-Arctic ice cod Arctogadus glacialis, which was found in almost half of the samples, is known to be most abundant at 300−400 m (Mecklenburg et al. 2018). The species is not common in Disco Bay (Jordan et al. 2003), so we regard it as a new discovery that A. glacialis is abundant in the Icefjord.

The absence of apex predators; favorable conditions for Greenland halibut and ringed seals
It is locally known that both Greenland halibut and ringed seals are generally larger inside the Icefjord compared to those caught in Disco Bay (Schiøtt et al. 2021). The reason for this difference might be due to (1) absence of the Greenland shark, (2) absence of polar bears, and (3) fishing and hunting activities being limited to the winter and spring. The Greenland shark is an abundant species around Greenland, often caught as bycatch in the fisheries, with highest densities in the Disco Bay area (Mac-Neil et al. 2012, Nielsen et al. 2014. Despite this, none of the informants reported the Greenland shark as bycatch in the Greenland halibut fisheries in the Icefjord (Table S1). Based on local knowledge and the lack of eDNA in the current study for this species, we can assume that the Greenland shark most likely does not reside in the Icefjord. The absence of this species could potentially be a contributing factor in explaining why the ringed seals and Greenland halibut that reside in the Icefjord are larger compared to those caught in Disco Bay, as the Greenland shark is considered an apex predator that also actively feeds on live seals and large Greenland halibut (Nielsen et al. 2014). A more likely explanation is that an intensive Greenland halibut fishery in Disco Bay is effectively removing all medium-sized fish, so they do not have the chance to grow large (Fredenslund 2022), where as in the Icefjord the less intensive fisheries make it more likely that the Greenland halibut can grow large before capture. Another apex predator that can be considered absent in the Icefjord is the polar bear, as several years can go by without sighting polar bears in the Icefjord (local knowledge). Fewer than 10 polar bears are caught each year in Disco Bay, while fewer than 3 bears are caught around Ilulissat (Grønlands Statistik 2016), and it is assumed that they are caught well before they reach the Icefjord. Predation from polar bears on ringed seals would have required at least 43 seals annually per bear to sustain their survival (Kingsley 1998) and their absence would relieve a large number of seals in the Icefjord from predation. The third limited predation factor in the Icefjord is from human activities, which are limited to the winter and spring (Schiøtt et al. 2021), giving both seals and fishes a refuge from human predation from late spring until late autumn inside the Icefjord. The limited seasonal (fishing/hunting) activities or total absence of apex predators (polar bear and Greenland shark) could suggest that the Icefjord is a refugium during part of the year for both ringed seals and Greenland hali but that reside in the Icefjord, placing them at the top of the local food chain, allowing them to grow older and larger.
Their larger size could also be the result of higher prey abundances, as more abundant prey positively influences predators' physical response through faster rate of growth, reproduction, and lower mortality (Asseburg et al. 2006). We currently have no information about the abundance of zooplankton and prey fish inside the Icefjord, but studies from other areas in Greenland have shown that plankton communities differ significantly between offshore areas and the inner fjords (Munk et al. 2003, Arendt et al. 2010 and are influenced by sea ice conditions (Møller & Nielsen 2020). Whether this is also the case in the Icefjord is an important question for future research in the area.

Polar cod vs. capelin; indications of different water types
When assuming that higher read count equals higher abundance, it can be assumed that polar cod is more abundant in the southern and northern branches compared to the main fjord, as polar cod is highly represented in the eDNA read counts in the samples collected both in the southern and northern branches of the Icefjord, while capelin eDNA have higher relative read counts -and are therefore most likely more abundant -in the main fjord (Fig. 4).
Local fishers and hunters know that harp seals forage on shoals of capelin in the ice-free areas of the main fjord (Table S1). This knowledge was supported by a stomach content analysis of harp seals caught in the Icefjord, which showed that 78% of the stomach content consisted of capelin (S. Schiøtt & A. Rosing-Asvid unpubl. data). This suggests that the main fjord is more influenced by inflows from Disco Bay (deep SPMW when inflows are warm and BBPW [see Fig. 2] when inflows are cold; Rysgaard et al. 2020) and therefore has higher abundance of capelin. The northern and southern branches have more coastal standing water masses (Fig. 2C), where polar cod is somewhat more abundant (Fig. 4). The distribution of capelin and polar cod in the different areas of the fjord system appears to follow the different water types seen in Fig. 2C, which resemble that shown in Hop & Gjøsaeter (2013), where capelin is associated with Atlantic water masses (Raymond & Hassel 2000, Rose 2005, while polar cod is associated with cold polar water. The TS relationship in Fig. 2C also suggests that the southern branch is different from the northern branch during the summer. The northern branch has CW masses deeper down the water column during July 2019, while the southern branch has completely different water types during the same period (Fig. 2C). The δ 18 O analysis also shows that the southern branch is slightly more influenced by meltwater from the ice sheet compared to the northern branch ( Fig. 3; y-intercept close to values of ice sheet δ 18 O), suggesting that the shallower sill here reduces the water exchange. This suggestion is supported by higher salinities below 40 m depth in the southern branch during July 2019 ob served in Fig. 2C, which might be the remnants from a previous dense inflow to the southern branch from the main fjord/Disco Bay that did not reach the northern branch. Unfortunately, no CTD measurements from around Ilulissat or from the Greenland Ecosystem Monitoring Program (www.g-e-m.dk) were made during the period April−July 2019 to support this hypothesis. TS analysis of the water types (Fig. 2C) suggests that the water masses in the Icefjord are a mixture of SPMW and CW. It seems that there was an inflow of water masses into the northern branch around February 2020, giving higher salinities in the deeper layers (Fig. 2C), while a different inflow also occurred in the southern branch sometime during the period April−July 2019, giving higher salinity below 25 m depth (Fig. 2C). In addition to these differences, it should be noted that there usually is an open water polynya at the entrance to the southern branch during winter and spring (local knowledge, supported by satellite imagery), which could have an influence on the different communities observed. These different conditions could contribute to the explanation of the observed higher diversity of taxa in the southern branch compared to the main and northern branch of the fjord system (Table 1, Fig. 5): taxa could be carried by the inflow of water masses into the southern branch affected by the Coriolis force that forces the inflow of water masses clockwise and therefore towards the southern branch first; alternatively, it may be the result of higher productive conditions caused by the open water polynya at the entrance to the southern branch. The lower number of taxa observed in the samples collected at 5 m depth in the northern branch during July 2019 (Table 1, Fig. 5) could be the result of lower salinities in the upper 40 m (Fig. 2B), which could have resulted in a displacement of the community of fish, forcing them to a greater depth where higher salinities are found. These lower salinities in the upper 50 m in the northern branch are considerably deeper than in the southern branch, where similarly low values are only observed in the upper 10 m of the water column (Fig. 2B). These could partly explain the observed higher diversity of taxa in the southern branch in July 2019, as the fish might not have been as strongly displaced by the lower salinities observed in the upper 10 m in the southern branch, resulting in the higher number of taxa observed in the samples from 5 m depth compared to the samples in the northern branch. The lower salinities and lower temperatures found around 50 m in the northern branch are likely the results of a greater presence of large icebergs that traveled northward from the main fjord. In the field, it was noted that there were no large icebergs in the southern branch during July 2019 compared to the northern branch, likely caused by the shallow sill at the entrance to the southern branch blocking the entrance of large icebergs from the main fjord. The δ 18 O values in Fig. 3 also indicate that the stratification of the water masses is strongest during the summer due to a higher influx of meltwater from the glaciers, resulting in lower δ 18 O values in the upper layers in the water column. During winter and spring, the water masses are well mixed, resulting in more clustered δ 18 O values (Fig. 3) and less stratification (Fig. 2), which supports the CTD measurements showing differences between locations and between seasons that could explain the observed differences in species abundance, seen in Fig. 6.

Marine mammals
Harp seals migrate to their breeding grounds around February (local knowledge, Rosing-Asvid 2010) and pass the Icefjord on their way down the west coast of Greenland in large numbers, supporting the much higher proportional read count contribution for this species in the samples collected in February 2020 (Fig. 4). Even though the ringed seals are year-round residents in the Icefjord (Schiøtt et al. 2021, A. Rosing-Asvid et al. unpubl. data), much lower read count proportions are observed for this species compared to harp seals. The difference can be explained in their foraging behaviors; adult harp seals forage in large herds, whereas ringed seals forage individually (Rosing-Asvid 2010), which could explain the relatively lower proportion of ringed seal eDNA compared to harp seals (Fig. 4). Local informants reported that whales are rarely sighted inside the Icefjord, which is also observed in their eDNA relative read counts that are much lower compared to ringed seals and harp seals (Fig. 4). Common minke whale, fin whale, and narwhal eDNA are observed in low relative read counts (Fig. 4) and only found in the samples from July 2019 (Table 1). These whales are common on the west coast of Greenland, and therefore it is no surprise to see their presence, but as their proportional contributions are somewhat low, their presence can be considered low, which supports the local knowledge that whales are rarely sighted inside the Icefjord (Table S1). The beluga, on the other hand, has high relative read counts in the samples collected in the main fjord (Fig. 4), so it is more likely that they were present in higher numbers in the main fjord during December 2019 (Fig. 4).

Ice conditions; future expectations
As mentioned above, the locals reported that the ice-clearing events in the main fjord have become more frequent and that the icebergs have become smaller. These changes might have cascading effects on the local ecosystem, as changes in the ocean temperatures, salinities, and ice conditions have a direct effect on lower trophic levels (Kortsch et al. 2012, Foss heim et al. 2015, Mueter et al. 2020, therein affecting local fish species (Hop & Gjøsaeter 2013, Bouchard et al. 2017, Bouchard & Fortier 2020, Møller & Nielsen 2020, Mueter et al. 2020. Thus, it is highly relevant to document changes in the local ecosystems, making eDNA studies an at tractive means of monitoring marine vertebrates in the fjord system. This study can be used as a starting point in doing so.

CONCLUSIONS
The aim of this study was to apply eDNA metabarcoding of seawater samples to enhance our understanding of the understudied ecosystem of Ilulissat Icefjord, focusing on spatial and seasonal changes in the species composition of marine vertebrates. We complemented the eDNA approach with local knowledge and hydrographic observations. The results from eDNA metabarcoding and hydrographic observations fit well with the knowledge gathered from local hunters and fishers on seasonal changes in the presence of e.g. Arctic char and harp seals as well as differences in the fish communities in the different areas of the fjord system. Results from eDNA metabarcoding alone would not have provided the same detailed description of the seasonal and spatial differences of these communities. This highlights the advantage of an interdisciplinary approach and collaboration with locals, especially in areas that are difficult to access. The Icefjord can be considered a refugium for both the ringed seals and Greenland halibut for part of the year, and a change in the fjord's accessibility would change their apex position, which highlights the need for monitoring changes in the fjord system. As trawl surveys or surveys from research vessels are currently not possible to conduct in the Icefjord due to heavy ice conditions, eDNA metabarcoding and hydrographical observations are currently the best options available to monitor the presence of fish and marine mammals inside the Icefjord and to look for changes in the water masses of the fjord system.