This study was conducted in the Bering Sea and the Chukchi Sea in the Arctic Ocean during a cruise (MR20-05C) conducted by the ‘R/V Mirai’ of the Japan Agency for Marine-Earth Science and Technology from September 19 to November 2, 2020 (Figure 1). The Bering Sea is surrounded by Alaska to the east, Siberia to the north-west, the Kamchatka Peninsula to the west, and Aleutian Islands to the south. The Bering Sea is divided into two regions: the southwestern part is an ocean basin with a depth of approximately 4,000 m, and the northeastern part is a continental shelf with a depth of less than 200 m. The Chukchi Sea, which connects to the Bering Sea via the Bering Strait, is composed of a continental shelf with a depth of approximately 50 m and is adjacent to the Canada Basin with a maximum depth of approximately 4,000 m.

Several distinct water masses have been found in the study area (Danielson et al., 2017). The Bering Slope Current, which flows northwestward along the continental slope, forms the boundary between the southwestern and northeastern parts of the Bering Sea. Three water masses have been recognized in the Bering Sea (Coachman et al., 1975): The Anadyr Water (AnW) is distinguished by its relatively saline, cold, and nutrient-rich water; The Alaska Coastal Water (ACW) is distinguished by its relatively warm and less saline water; The Bering Shelf Water (BeSW) having intermediate characteristics between AnW and ACW. These water masses flow into the Chukchi Sea through the Bering Strait and comprise the Bering-Chukchi Waters (BCSW; Bering-Chukchi Summer Water and BCWW; Bering-Chukchi Winter Water; Danielson et al., 2017). After entering the Chukchi Sea, these water masses flow northward along three major conduits: the Herald Canyon in the west, the Barrow Canyon in the east, and the Central Channel across the mid-shelf (Danielson et al., 2017). The Beaufort Gyre circumscribes the northern shelf slope of the Chukchi Sea shelf in a clockwise direction. The surface of the Canada Basin is covered by a polar mixed layer consisting of a low-salinity sea ice meltwater layer (MW) with near-freezing temperature (Danielson et al., 2017; Skjoldal, 2022).

The cruise was conducted during the open-water season. According to satellite-based sea ice concentration (SIC) data from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data, Version 2, with a 25 km grid (provided by the National Snow and Ice Data Center (NSIDC); https://nsidc.org/data/data-access-tool/NSIDC-0051/versions/2), sea ice covered the Alaskan side of the Canada Basin. The northernmost edge of our study area, located around the Chukchi Borderland, was in the marginal ice zone (Figure 1).

The SIC data was also used to approximate the distance between the sampling sites and the sea ice edge. The daily SIC data corresponding to the study period were downloaded from the previously mentioned NSIDC website. The data were visualized using panoply version 5.2.3 (available online: https://www.giss.nasa.gov/tools/panoply/) and then exported as a KMZ file. These KMZ files were imported to Google Earth Pro 7.3.6.9345 (available online: https://www.google.co.jp/earth/about/versions/). The sea ice edge was defined as the contour of SIC = 0.15. The shortest distance between the sampling sites and the sea ice edge of the day was then calculated using the “Measure” function in Google Earth Pro. To avoid the path on a map intersecting land, the coordinates of Diomede, located in the middle of the Bering Strait (65° 45’ 30” N, 168° 57’ 06” W), and a sampling site located off Point Hope (68° 30’ 1.31” N, 168° 49’ 41.84” W), were used as waypoints when calculating the distance from the sampling sites in the Bering Sea.

eDNA samples were collected from 29 sites in the Chukchi Sea and nine sites in the Bering Sea from October 3 to October 24 (Figure 1). The geodesic distance from every sampling site to the nearest one ranged from 2 to 255 km (median [first quartile, third quartile] = 77 [24, 149] km). Based on the distance between our sampling sites and the current speed (m/s) that was obtained with a hull-mounted Acoustic Doppler Current Profiler (Teledyne RD Instruments, San Diego, CA, USA), it is estimated that eDNA shed at a particular site would take approximately 0.1–32 days (3.8 [1.6, 12.9] days) to drift from one sampling site to the nearest neighboring site.

As the surface water, seawater pumped from approximately 4.5 m depth was collected from the laboratory faucet of the ship. Because the seawater was used for continuous measurements of environmental conditions and water quality, it was filtered using a coarse metal scrub brush along the path of intake to remove large materials that could interfere with the measurements. We maintained a constant flow (approximately 5 L/min) of seawater throughout the cruise to ensure that the seawater in derivation pipes and tubes was always replaced with newly taken seawater. The insides of the pipes and tubes were considered to be constantly rinsed with fresh seawater. In addition to the surface water, water samples were collected from the middle or bottom layers using 12 L Niskin bottles mounted on a conductivity–temperature–depth (CTD)/Carousel Water Sampling System at six stations (Figure 1). The middle layer was defined as beneath the mixed layer or at the depth with the highest light attenuation rate, which is expected to be the most turbid layer (Supplementary Figure 1). The bottom layer was defined as being 5 m above the sea floor. For each water sample, a clean, foldable plastic bag was filled with 12 L of collected water.

Detailed sampling procedure for eDNA was described in Kawakami et al. (2023). Briefly, water samples were vacuum-filtered onto a Sterivex-HV filter unit (0.45 µm pore size, PVDF; Millipore, Burlington, MA, USA) using an aspirator (GAS-1N; AS ONE, Osaka, Japan) connected to a manifold (Multivac 310-MS; Rocker, Kaohsiung City, Taiwan) in the ship’s laboratory, which was not used for processing organism samples. Filtration was halted when the filter became severely clogged. The filtration volume was calculated by subtracting the volume of unfiltered water from the initial volume (12L). Consequently, the filtration volume ranged from 1.6 to 12 L. The filtered sample was immersed in 1.6 mL of RNAlater (Thermo Fisher Scientific, Waltham, MA, USA) and stored at –20°C until DNA extraction. For every 10–12 samples of seawater, negative controls were created by filtering 500 mL of Milli-Q water following the same procedure as seawater sampling.

To prevent cross-contaminations, all equipment used for water collection and filtration was decontaminated by soaking them in a 2% bleach solution (approximately 0.12% NaOCl) for at least 10 min. To remove any remaining bleach, the equipment was subsequently rinsed with tap and Milli-Q water. The workspace was decontaminated using bleach disinfectant spray (approximately 0.1% NaOCl). For each water collection and filtration, new gloves and masks were used.

Temperature (°C), salinity (psu), and chlorophyll a (chl a) fluorescence (RFU) of the surface water at the sampling sites were obtained from a Continuous Sea Surface Water Monitoring System (Marine Works Japan Co., Ltd., Kanagawa, Japan). Their vertical profiles were measured by the CTD at sampling sites where the middle and bottom waters were collected. The CTD data are available in the Arctic Data archive System (ADS) managed by the National Institute of Polar Research, Japan, under the ID of A20220803-010 (Nishino, 2020; Nishino, 2022).

DNA was extracted from each filter using the Qiagen DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany), according to Wong et al. (2020), who modified the original extraction protocol to improve the extraction efficiency from the Sterivex filter. Extraction was performed using vacuum-driven liquid processing with a manifold (QiaVac; Qiagen) and vacuum pump (MAS-01; AS ONE). Throughout the process, a VacValve and VacConnector (Qiagen) were utilized in conjunction with the manifold to reduce the risk of contamination from backflow into the columns and direct contact with the manifold.

After removing the RNAlater, 2 mL of a lysis buffer mix containing PBS (990 L), Buffer AL (910 L), and Proteinase K (100 L) was introduced to the Sterivex through its outlet using a 2.5 mL syringe (Terumo Corporation, Tokyo, Japan). The filter was sealed with luer fittings and incubated at 56°C without rotation for 30 minutes. Subsequently, the filter was centrifuged at 2500 g for 10 minutes, and the lysate was retrieved in a 3 mL polypropylene tube. The lysate was filtered using a DNeasy Blood and Tissue kit column after being mixed with 1 mL of 99.5% ethanol (FUJIFILM Wako Pure Chemical Corporation, Osaka, Japan). The column was then washed sequentially with 0.8 mL of AW1 buffer and 0.8 mL of AW2 buffer, followed by 2 minutes of centrifugation at 21,300 g for drying. The adsorbed DNA was eluted from the column into a DNA LoBind Tube (Eppendorf, Hamburg, Germany) twice with 75 μL AE buffer to maximize recovery. A total of 150 μL DNA extract was collected from each filter and stored at –20°C.

We designed novel primers and probes that exclusively amplify DNA derived from polar cod in samples collected from the Bering and Chukchi Seas. The mitochondrial D-loop region was chosen as a candidate region for designing a new specific marker because it is one of the major target regions adopted for species-specific eDNA detection owing to its high inter-specific variability (Tsuji et al., 2019).

Twelve gadid species are known to be distributed in the Arctic region (Mecklenburg et al., 2018). We searched and downloaded the D-loop sequences available for these species from GenBank (http://www.ncbi.nlm.nih.gov/), resulting in a set of D-loop sequences for polar cod and eight other species (saithe, Pollachius virens; Arctic cod, Arctogadus glacialis; saffron cod, Eleginus gracilis; Atlantic cod, Gadus morhua; walleye pollack, G. chalcogrammus; Pacific cod, G. macrocephalus; Haddock, Melanogrammus aeglefinus; and whiting, Merlangius merlangus).

Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used to generate specific primer pairs for polar cod on the target region. For running Primer-BLAST, a sequence of the D-loop region of polar cod extracted from the reference sequence (NC_010121) was provided as a PCR template, and those of non-target species were set as a database for checking primer pair specificity. T_m of the candidate primers was calculated using the T_m calculator provided by ThermoFisher Scientific (https://www.thermofisher.com/jp/ja/home/brands/hermos-scientific/molecular-biology/molecular-biology-learning-center/molecular-biology-resource-library/hermos-scientific-web-tools/tm-calculator.html). Subsequently, primer pairs with an amplicon length of less than 180 bp and a T_m of 58-60°C were selected. Following that, an MGB probe was designed on the amplicon sequences using PrimerExpress 3.0 (Thermo Fisher Scientific, Waltham, MA, USA) using the default settings.

The specificity of the newly designed primers and probe for polar cod was confirmed in silico using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) against the nr database. Consequently, we selected the primer set and probe described in Table 1 and Supplementary Figure 2 for the polar cod eDNA assay. All available sequences for these regions of polar cod were identical, except for one sequence containing one mismatch in the probe region.

The specificity of the assay was further confirmed in vitro using DNA samples extracted from the tissues of five gadid species (polar cod, Pacific cod, walleye pollock, Atlantic cod, and saffron cod), which may have co-occurred in the study area. Their DNA was extracted using the DNeasy Blood & Tissue kit (Qiagen) following manufacturer’s protocol from the following tissues: the ethanol-preserved muscle of polar cod collected in the Chukchi Sea during the cruise of T/S Oshoro-Maru of Hokkaido University; the ethanol-preserved testis of pacific cod captured off Hokkaido, Japan; purchased dried fish for walleye pollock and saffron cod captured off Hokkaido, Japan; a dried dorsal fin of Atlantic cod, collected around Norway. In the specificity test, these extractions were diluted 10 pg/µL and 100 pg/µL after measuring the concentration using a Qubit dsDNA HS Assay Kit and Qubit 4 Fluorometer (ThermoFisher Scientific).

qPCR with a newly designed species-specific eDNA assay was performed in a StepOne plus real-time PCR system (Life Technologies) with a 15 µL reaction containing 1x Environmental Master Mix 2.0 (ThermoFisher Scientific), 900 nM of each primer, 125 nM of a probe, 0.075 μL of AmpErase™ Uracil N-Glycosylase (Thermo Fisher Scientific, Waltham, MA, USA), and 2 μL of DNA template. A no-template control (NTC) was included in the PCR run using 2.0 μL of nuclease-free water (UltraPure; Invitrogen, Waltham, MA, USA) in place of the DNA template. The thermal conditions were as follows: 50°C for 2 min, 95°C for 10 min, 55 cycles at 95°C for 15 s, and 60°C for 1 min.

A newly designed assay was applied to eDNA samples collected from the Bering and Chukchi Seas to estimate the distribution of polar cod. The assay was performed under the same PCR conditions as those described above. A 2 µL of eDNA extract was used as a template in each PCR. To evaluate cross-contamination, a no-template control (NTC) was included in each PCR run using 2 μL of nuclease-free water (UltraPure; Invitrogen, Waltham, MA, USA) instead of the DNA template. Three technical replicates were used for each sample.

To determine the quantity of target DNA in each sample, 10-fold serial dilutions ranging from 30 to 30,000 copies/reaction of the artificially-synthesized double-stranded DNA (gBlocks; Integrated DNA Technologies, Inc., Coralville, IA, USA) containing the target sequence were included in triplicate for each qPCR run. The concentration of the stock standard solution was quantified using a Qubit dsDNA HS Assay Kit and a Qubit 4 Fluorometer (Thermo Fisher Scientific) and then diluted in TE buffer. To eliminate suspicious detection, reaction resulting in < 1 copies/reaction was converted to 0. After conversion, the mean of the triplicates was used to represent the number of eDNA copies in each sample. The eDNA concentration in seawater (copies/L) was calculated by converting the number of eDNA copies per reaction to the total number of eDNA copies present in the DNA extract, and then dividing this value by the filtration volume.

For DNA extraction and PCR preparation, filter samples were processed in a laboratory dedicated to eDNA analysis, which was isolated from a room for PCR, ensuring neither tissue nor PCR products were handled. PCR amplification was performed in a separate room. During every step of the molecular work, researchers always wore new gloves and masks. The laboratory and all equipment (micropipettes, incubator, and centrifuge) were cleaned with DNA-OFF (Takara, Shiga, Japan) and distilled water. In addition, the extraction apparatus (manifold, valves, and tubes) was decontaminated with a bleach solution, as detailed previously.

All statistical analyses were performed using R v.4.1.2 (R Core Team, 2021) in combination with RStudio v.2022.02.3 (available online: https://www.rstudio.com/) and the packages implemented in it. Additionally, GMT 5.4.5 (Wessel et al., 2013) was used to visualize the results.

To explore the association between the study sites and water masses, the sampling sites were classified using clustering and ordination analyses based on water temperature, salinity, and chl a, which are potential indices of distinct water masses across the Bering and Chukchi Seas (Danielson et al., 2017). First, we performed a clustering analysis based on Euclidean distances using Ward’s linkage. The similarity profile (SIMPROF) test (Clarke et al., 2008) using the “simprof” function of the “clustsig” package in R was subsequently applied to identify statistically significant clusters (the significance threshold was set at 0.05). The SIMPROF test examines the multivariate structure within groups of samples by comparing the similarity profile constructed from the results of the cluster analysis with that constructed from the original data using the iterative permutation procedure. The results of the cluster analysis were visualized by Principal Component Analysis (PCA) using the “prcomp” function of the “stats” package and “fviz_pca_biplot” function of the “factoextra” package in R. Chl a was included in these analyses after log transformation. Before these analyses were conducted, the environmental parameters were normalized.

The environmental conditions were compared between the sampling sites where polar cod eDNA was detected or not detected using the Wilcoxon rank-sum exact test after assuming homoscedasticity based on the results of the Bartlett test. The significance level was set at 0.05.

Multiple logistic regression analysis was used to relate the detection of polar cod eDNA to environmental conditions (water temperature, salinity, chl a, and topography). Prior to the logistic regression, the qPCR results were converted to presence/absence data. The distance from the sea ice edge was excluded from the analysis due to a significant and strong correlation with water temperature (Pearson’s product moment correlation: r = 0.906, t = 12.823, p < 0.001). Chl a was incorporated into the model after log transformation. The topography was roughly determined based on the bottom depth as shelf (< 200 m) and slope or basin (≥ 200 m). The logistic model was fitted to the data with a binomial error distribution and a logit link using the “glm” function in the “stats” package in R. After constructing a full model, the best model was selected using the “stepAIC” function in the “MASS” package in R based on Akaike’s information criterion (AIC).

As a result of qPCR using DNA extracted from the tissues of five gadid species, including polar cod, our newly designed assay exclusively detected the DNA of polar cod at both concentrations (10 pg/µL and 100 pg/µL), except for one sample of walleye pollock containing 100 pg/µL of DNA (Supplementary Figure 3). Although nonspecific detection of walleye pollock DNA was observed, the DNA was not detected at lower concentration (10 pg/µL) and its Ct (45.98) was clearly higher than the intercept of standard curve (41.03, Supplementary Table 1), which corresponds to 1 copy/reaction. This suggested that the risk of false detection of walleye pollock can be negligible when using 1 copy/reaction as a cut-off value for quantity. Therefore, this assay was assumed to be sufficiently specific and was used in subsequent analyses.

The R², intercept, log-linear slope, and efficiency of the standard curve across four runs for eDNA samples were 0.998 ± 0.001, 41.03 ± 0.37, -3.50 ± 0.07, and 93.19 ± 2.47 (mean ± SD), respectively (Supplementary Table 1). None of the negative controls (field negative controls and NTC for PCR) yielded positive results for polar cod eDNA, indicating that cross-contamination during the experiment was negligible and that the qPCR results obtained in this study were reliable.

The sampling sites were classified into five significant clusters based on the SIMPROF results (Supplementary Figure 4). Overlaying the clustering results on the PCA ordination and a geographical map revealed that these clusters corresponded well with the known water masses throughout the study area (Figure 2A, Table 2, Supplementary Figure 5). Cluster 1 consisted of sampling sites south of 62°N in the Bering Sea, where the water was saline, warm, and high in chl a concentration. The sampling sites across the Bering Strait and central Chukchi Sea shelf comprised two distinct clusters (clusters 2 and 3). Although salinity was comparable between these three clusters, both temperature and salinity decreased gradually with increasing latitude (Table 2, Supplementary Figure 5). The sites along the northernmost edge of our study area (approximately 75°N) were also classified as having subzero temperatures and low salinity (cluster 5). The sites near the Point Barrow were classified as another cluster (cluster 4) with intermediate temperature and salinity compared to those of clusters 3 and 5. In addition, both water temperature and salinity tended to decrease with increasing proximity to the sea ice edge. These clusters corresponded to the water masses identified by Danielson et al. (2017) based on temperature and salinity: ACW (cluster 1), BCSW (clusters 2 and 3), and MW (clusters 4 and 5), respectively.

Polar cod eDNA was detected primarily in the central part of the Chukchi Sea shelf and in the northernmost observation line, which was located on the shelf slope off the Point Barrow and marginal ice zone (Figure 2B). In contrast, polar cod eDNA was rarely found in the Bering Sea. Based on the classification of water masses, polar cod eDNA was most frequently and abundantly detected in the sites classified as clusters 4 and 5, which showed cold, less saline water (92% with 46 ± 69 copies/L in cluster 4, and 100% with 156 ± 289 copies/L in cluster 5; Table 2). In contrast, no polar cod eDNA was detected in the sites belonging to cluster 1, which showed warm, saline water. At sites belonging to clusters 2 and 3, polar cod eDNA was detected sporadically at moderate concentrations (50% with 7 ± 11 copies/L in cluster 4, and 40% with 95 ± 193 copies/L in cluster 5).

Polar cod eDNA was also detected in the middle and bottom layers of all six sites examined for vertical comparison, including sites where no polar cod eDNA was detected in surface samples (Figure 3). Despite the different environmental conditions of the site’s surface water, polar cod eDNA was detected in the bottom layer across the Chukchi Sea shelf at depths of 50–120 m.

The juxtaposition of the concentration of polar cod eDNA on the temperature–salinity diagram showed that the distribution of polar cod eDNA in the surface water was clearly limited to < 5°C in temperature and < 32 psu in salinity (Figure 4). Maximum concentrations of eDNA in surface water occurred at approximately 0°C and 30 psu. There were significant differences in water temperature and salinity between sites where polar cod eDNA was detected and those where it was not (Supplementary Figure 6; Wilcoxon rank-sum exact test; W = 35, p < 0.001 for water temperature; W = 32, p < 0.001 for salinity). Depending on the depth, the relationship between eDNA distribution and environmental conditions varied; polar cod eDNA was detected in more saline water in the middle or bottom layers (> 32 psu) than in the surface water.

After developing a logistic regression model relating the presence or absence of polar cod eDNA to four environmental conditions (water temperature, salinity, log-transformed chl a, and topography), the model containing only water temperature as an independent variable was deemed the best (Supplementary Table 2). For each 1°C increase, the odds of detecting polar cod eDNA decreased by 50% (95% CI [24%, 67%]). The probability of detecting polar cod eDNA decreased with increasing water temperatures (Figure 5). The water temperature with a 50% detection probability was predicted to be 3.7°C.

As eDNA analysis in marine ecosystems detects DNA in water, which is a trace of individual organisms, it can only indirectly suggest the presence of target species. Various biotic factors, such as species, physiological conditions, body size, and developmental stage, influence the amount and rate of eDNA shedding (e.g., Maruyama et al., 2014; Klymus et al., 2015; Takeuchi et al., 2019). The ecology of eDNA following its release is also affected by degradation, particle sinking, adsorption to sediments, and dispersion (Barnes and Turner, 2016; Stoeckle et al., 2017). False negatives (species not detected where they are present) and false positives (target species absent but DNA recovered) in eDNA detection can also arise from a variety of biological, physical, and chemical factors, and technical issues (Barnes and Turner, 2016; Cristescu and Hebert, 2018). Consequently, it is often difficult to infer the distribution and abundance of the source organism based on the detection and quantification of eDNA. These obstacles could confound our inferences about the distribution and habitat preferences of polar cod in the present study, necessitating a careful inspection of the results.

A properly designed primer pair and probe are critical to avoid false results. Primers and probes that lack sufficient specificity can amplify unintended species, whereas high intraspecific variation in priming sites can result in failure to recover eDNA from the target species. Therefore, we rigorously tested the performance of the new assay to ensure its specificity for polar cod, employing publicly available sequences of polar cod and other related species, as well as PCR using DNA extracts of gadoids that may have co-occurred with polar cod in the study area. The results of in silico and in vitro specificity tests confirmed that the newly designed assay had the desired specificity, i.e., the assay was expected to amplify polar cod DNA exclusively. However, the risk of false results remains due to the uncertainty regarding unknown intraspecific diversity. A recent study on the genetic structure of polar cod populations using microsatellite loci found that polar cod in the North Bering, Chukchi, and Western Beaufort Seas formed a genetically distinct group from other regions, with very little genetic differentiation within the group (Nelson et al., 2020). A study that used mitochondrial DNA also suggested no significant structure from the northern Bering Sea to the western Beaufort Sea (Wilson et al., 2020). Although these results cannot directly imply negligible intraspecific diversity in our target region (mitochondrial D-loop), it is unlikely that the new assay produced false negative results due to intraspecific diversity.

As fish eDNA in the ocean is thought to represent only a small fraction of the DNA in seawater (Stat et al., 2017), the concentration of eDNA may be influenced by the dispersion and degradation processes critical for its detection. Therefore, insufficient sampling effort can lead to false-negative errors. However, the volume of water filtered in this study (mean ± SD, 7.4 ± 3.2 L per filter) exceeded that employed in previous fish eDNA surveys in the ocean to comprehensively detect the local diversity of fish fauna (0.5–2 L of seawater per sample, e.g., Thomsen et al., 2016; Sigsgaard et al., 2017; Afzali et al., 2021). Given the prevalence of polar cod in the Arctic and a large amount of water filtered in this study, we expected the occurrence of false negative errors to be negligible.

False positives caused by contamination of space and time should also be considered. Although dispersion lowers eDNA concentration, it can increase the probability of detecting non-local and non-temporal eDNA originating from spatially and temporally remote sources from a study site (Hansen et al., 2018). In the Antarctic, eDNA from the ocellated icefish, Chionodraco rastrospinosus, is expected to disperse up to 648 km from its source in a detectable state, assuming a current velocity and slow decay rate in subzero water temperature (Cowart et al., 2018). Bottom sediments and sea ice are other potential eDNA sources in the ocean, transporting allochthonous extracellular DNA (Collins and Deming, 2011; Kuwae et al., 2020). However, these cases likely represent extremes. Half-lives of marine fish eDNA have been estimated between approximately 7 and 60 h under 4–40°C and 71 h under -1°C (Collins et al., 2018). A field experiment in the temperate zone showed that eDNA shed by white trevally (Pseudocaranx dentex) was only detectable within a narrow range (<100 m) in a bay with a weak current (Murakami et al., 2019). These studies implied that eDNA detection is thought to represent a fish community on a short temporal and small spatial scale. In the present study, approximately 3.8 days (91 h) in the median were required for eDNA shed from a particular site to be transported to the nearest site, assuming the passive transportation of eDNA by the current. Although there is a possibility that the degradation of eDNA could be slower than expected under the Arctic environment, it could be presumed that our sampling sites did not share the same source for eDNA due to the temporal distance between sites.

The distribution of polar cod eDNA corresponded well with the putative water mass distribution, suggesting that the polar cod was heterogeneously distributed in the study area and was likely affected by site-specific environmental conditions. Although eDNA may reliably suggest the polar cod distribution as discussed above, when considering interactions with various environmental conditions, the complex relationship between eDNA concentrations, abundance, and body size of the source organism must be clarified to accurately convert the eDNA distribution into reliable data for conservation purposes (Collins et al., 2018; Jo et al., 2019; Jo and Minamoto, 2021). Across various taxonomic groups, rearing experiments have generally shown positive correlations between the amount of eDNA and the abundance, body weight, and developmental stage of source organisms (e.g., Takahara et al., 2012; Minegishi et al., 2019; Takeuchi et al., 2019). This relationship between eDNA concentration and the abundance of target species has also been supported in coastal environments (Thomsen et al., 2016; Yamamoto et al., 2016; Sato et al., 2021), with some exceptions in the open ocean (Knudsen et al., 2019; Shelton et al., 2022). Although uncertainty remains regarding obtaining comparable quantification from the eDNA results, we discuss what and how environmental conditions shape the polar cod distribution, carefully considering that limitation.

Although oceanographic conditions vary temporally and spatially, previous studies that intensively described polar cod distribution in the Chukchi Sea during the open-water period have shown their prevalence distribution across the shelf, but it was heterogeneous and linked to characteristics of the water mass. Polar cod occupy a remarkably wide range of habitats; they were collected from all parts of the water column, from the coast to offshore (Goddard et al., 2014; Logerwell et al., 2015). Gradual expansion in body size range, in which the maximum body length increased from 66 mm to 260 mm (corresponding to age-0/1 and age-3+, respectively) from the beach to the shelf surface, nearshore bottom, and benthic shelf, suggested a habitat shift linked to growth (Logerwell et al., 2015). Polar cod is the most abundant fish species in the region, except on the shelf surface, where capelin (Mallotus villosus) and saffron cod tend to represent a large proportion of the fish stock (Norcross et al., 2010; Eisner et al., 2013; Goddard et al., 2014; Logerwell et al., 2015). Studies focusing on the early life stage distribution of polar cod suggested that their larvae and juveniles are abundant in cold water ranging between approximately -1.8–5°C (Norcross et al., 2010; Kono et al., 2016; Deary et al., 2021), while their occurrence in warmer water (2.6–11.2°C) was also reported (Eisner et al., 2013).

Hotspots of polar cod eDNA were observed in the central part of the Chukchi Sea shelf, off the Point Barrow, and in proximity to the Chukchi Borderland, while they were scarce in the northern Bering Sea, which corresponds to the southern limit of its distribution (Mecklenburg et al., 2018). These results are congruent with the synthesized distribution of the historical data described above and the recent acoustic-trawl survey targeting polar cod (De Robertis et al., 2017). Therefore, eDNA detection successfully inferred the distribution of the polar cod. In this study, we mainly collected surface seawater for eDNA detection, but eDNA was also detected in the water in the bottom or middle layers. A recent study showed that the distribution of eDNA is vertically different and can reflect the vertical distribution of marine fishes (Canals et al., 2021). Although the comparison between the different depth layers in this study is preliminary, eDNA can be expected to indicate the vertical distribution of polar cod.

The key finding of this study is that eDNA showed a clear linkage between the presence of polar cod and water temperature, implying that polar cod prefer cold water masses as a habitat. Polar cod eDNA was only found in water masses at lower temperatures (<5°C), and detection probability decreased with increasing temperature. This finding is consistent with the previously reported temperature range in which polar cod larvae and juveniles are abundant (Norcross et al., 2010; Kono et al., 2016; Deary et al., 2021). Thus, eDNA provides new empirical evidence that supports the general assumption that polar cod prefer cold water. Rearing experiments have also suggested that polar cod are most active and capable of growing at temperatures below 10°C (Laurel et al., 2016). However, salinity may not have been a decisive factor in polar cod habitat selection because we detected their eDNA in a wide salinity range (25–34 psu) in our study area. In addition to temperature and salinity, other environmental conditions, interactions with other species (predators, prey, and competitors), and internal factors (e.g., density, age, and fish condition) can influence fish distribution (Planque et al., 2011). Considering that polar cod inhabit abundantly in the cracks and crevices beneath the sea ice (Gradinger and Bluhm, 2004), the presence of sea ice itself could also influence the distribution of polar cod and, consequently, the probability of eDNA detection. However, due to the close relationship between the distribution of cold water masses and the presence of sea ice in our study area, it was difficult to determine whether the habitat preference of polar cod was for the water temperature or the habitat availability under sea ice. Although the habitat preferences of polar cod are extensively reviewed (Geoffroy et al., 2023), their details are still being discussed and warrant further investigation.

According to previous studies, the eDNA hotspot in the central Chukchi Sea shelf observed in this study could have captured an aggregation of larval or juvenile polar cods. It has generally been assumed that polar cod spawn in the winter (November to April) beneath the sea ice (Thorsteinson and Love, 2016). Pelagic eggs and larvae that hatch from spring to summer disperse over a broader geographical area (Aune et al., 2021; Deary et al., 2021). Assuming the Kotzebue Sound is a potential source of polar cod, as this area is unknown as a spawning ground in the Chukchi Sea, individual−based models in combination with particle tracking predicted that simulated polar cod hatching in the Kotzebue Sound was transported offshore to the north and to the west, concentrated in the central Chukchi Sea around 70°N (Deary et al., 2021). This estimated dispersal process was consistent with the eDNA distribution. However, the inability of eDNA to distinguish the biological parameters of the source organism (e.g., age, body size, and physiological condition) precludes further estimation of the factors controlling its observed distribution. Another limitation of eDNA is that it does not provide immediate information about target species, necessitating laboratory work after a field investigation. The combined use of the eDNA technique, and the reconstruction of life history traits using otolith microchemistry or tissue stable isotopes, and hydroacoustic surveys with high spatial coverage would expand our understanding of polar cod ecology considerably.

The application of eDNA techniques would likely benefit various studies on tracking the consequences of environmental change, ecological management and conservation, and policy-making decisions (Bohmann et al., 2014). This study demonstrated that species-specific detection of eDNA is an effective method to describe the distribution of polar cod, providing baseline data to reveal its environmental impact. However, uncertainties from eDNA must be overcome before its broad adoption (Bohmann et al., 2014; Cristescu and Hebert, 2018). For example, we were unable to conclusively generate informative relative abundance and population structure data, due to the unknown relationship between the abundance or biological condition of the source organisms and eDNA detection, as well as the unknown ecology of eDNA in the Arctic. Currently, the abundance data of polar cod in the Chukchi Sea are limited (De Robertis et al., 2017; Geoffroy et al., 2023), and if the eDNA technique can indicate its abundance and population structure, it could provide a basis for projecting shifts in distribution and abundance. Experimental studies regarding eDNA shedding from polar cod and its persistence in a subzero environment would facilitate future applications of the eDNA technique. Furthermore, due to the inability of eDNA to provide biological parameters of the source organism, which is invaluable information for conservation and management purposes, the eDNA method should be considered as a complementary tool for other individual-based methods, such as target fishing, otolith microchemistry (Andrade et al., 2020), stable isotopes (Kohlbach et al., 2017), and life history analysis (Bouchard et al., 2017), rather than replacing them. The eDNA method could also supplement hydroacoustic surveys and observations using underwater optical camera by reducing ambiguity in species identification and enabling the production of distributional data at a large spatiotemporal scale with high taxonomic resolution.

Integrating eDNA techniques and oceanographic surveys will enhance our ability to detect predicted ecological changes in polar cod populations. eDNA samples can be collected concurrently with samples for biogeochemical, planktonic, and microbial studies, providing an additional layer for investigating biological interactions. A numerical hydrodynamic model incorporating oceanographical conditions (e.g., current, density, and temperature) and eDNA data (e.g., concentration and persistence) is a possible method for creating an accurate and quantified map of polar cod populations (Fukaya et al., 2020). In addition, recent advances in the application of eDNA, such as population genetic analysis (Sigsgaard et al., 2020; Tsuji et al., 2020) and analysis of eRNA associated with physiological conditions (Tsuri et al., 2021), may provide new insights into the population dynamics of polar cod. As this study may expand the applicability of eDNA techniques for Arctic studies, investigations targeting other species and accumulating basic knowledge of eDNA in the Arctic will facilitate the understanding of the processes and consequences of distributional shifts of fishes (Hollowed et al., 2013; Fossheim et al., 2015) predicted to occur under global warming in the future.