The overwintering ecology of fish in the seasonal ice zones (SIZs) remains largely unexplored owing to methodological limitations. Environmental DNA (eDNA) can reveal the distribution and diversity of fish species in various aquatic environments, thereby offering a possible solution to the methodological limitations of SIZ studies. Therefore, we aimed to detect the overwintering fish community in the ice-covered Saroma-ko Lagoon, located on the Okhotsk Sea coast of Hokkaido, and its inflow, using eDNA metabarcoding. eDNA extracted from under-ice seawater collected from the lagoon yielded 28 fish taxa, predominantly Clupea pallasii based on the relative DNA read abundance. Dissimilarity analysis suggested short-term temporal variations in eDNA composition in under-ice seawater, even at the same site. Nineteen fish taxa, predominantly Tribolodon brandtii and T. hakonensis, were detected in the eDNA extracted from under-ice river water. The high dissimilarity between eDNA results from under-ice seawater and river water suggested segregation of the overwintering community between the lagoon and river. Fish eDNA detected in meltwater from sea ice was assigned to five taxa, suggesting the entrainment of particulate matter containing fish eDNA during ice growth. The true species richness estimated based on eDNA results and discrepancies with historical reports suggested that sampling efforts need to be optimized for ice-covered environments to promote more comprehensive species detection. This study demonstrated the usefulness of using fish eDNA metabarcoding to study the ecology of overwintering fish under ice. Clarifying eDNA shedding patterns, persistence, and dispersal in under-ice environments would improve the reliability of this technique and expand its use in SIZs.

Study site

The Saroma-ko Lagoon is an enclosed coastal lagoon with a surface area of 151.59 km² and a maximum depth of 22 m. It is connected to the southern part of the Sea of Okhotsk by two inlets (Figure 1). The water mass of the eastern part of the lagoon consists primarily of Okhotsk Sea water, with freshwater input from the Saromabetsu River (Nomura et al., 2009). The sediments in the lagoon consist primarily of silt and sand (Nishihama & Hoshikawa, 1992). Almost the entire surface of the lagoon is covered with sea ice from early January to early April (Shirasawa et al., 2005). A weak but steady horizontal and vertical current (<0.02 m/s) likely driven by tidal fluctuations is observed beneath the sea ice (Nomura et al., 2024). The lagoon has moderate fish diversity (49 species) and is a composite of coastal marine fish (e.g., Clupea pallasii, Sebastes spp., Pleuronectidae spp.) and brackish water fish (e.g., Mugil cephalus) (Ministry of the Environment, Japan, 1993). Fish fauna include species targeted for under-ice fishing, such as Eleginus gracialis and Hypomesus nipponensis.

eDNA sampling

Sampling methods for eDNA collection are detailed by Nomura et al. (2022). Briefly, under-ice water was collected at the eastern part of the lagoon (St. L, which corresponds to the “eDNA site” in Nomura et al., 2022) and river mouth of the Saromabetsu River (St. R) on March 5 and 6, 2021 (Figure 1, Table 1). The water depth was approximately 12 m at St. L and 2 m at St. R. The distance between these two stations was approximately 4.25 km. The ice thickness at both sites was approximately 0.5 m. Water temperature and salinity in the lagoon ranged from -1.2 to -0.8 °C and 25.0 to 32.0 at St. L (Nomura et al., 2022), and -0.1 °C and 0.3 at St. R (Table 1), respectively. Surface sea water temperature off the Saroma-ko Lagoon in early March was -1.74 °C on average, with partial coverage of drifting sea ice (data was available in Japan Meteorological Agency: https://ds.data.jma.go.jp/kaiyou/data/db/kaikyo/series/engan/engandata_SP.html and Ice Information Center, Japan: https://www1.kaiho.mlit.go.jp/KAN1/drift_ice/ice_chart/calendar/2021/03.html).

Under-ice water was collected from a depth of 1.5 m (approximately 1 m below sea ice) using a water sampler or pump through an ice hole cut using an ice saw as described below. At St. L, water collection was repeated four times during the ebb tide (from 9:00 to 14:00) at approximately 1–2 hour intervals, considering short-term temporal variation in eDNA signals (Jensen et al., 2022). To evaluate the collection efficiency of under-ice seawater, two types of water samplers (a 1 L Kitahara sampler [Rigosha, Japan] and a 2 L Van Dorn sampler [Rigosha, Japan]) and two types of pumps (a hand pump [GH-0400D, Bosworth Company, East Providence, RI, USA] and an electric pump [YS-20, Sendak, Tokyo, Japan]) were used. At St. R, a peristaltic pump (Geopump, Geotech Environmental Equipment Inc., Denver, CO, US) was used to collect under-ice river water. Each collection filled at least two sterilized 10 L foldable plastic bags, and a benzalkonium chloride (BAC) solution (final concentration of 0.01%) was added immediately to prevent eDNA degradation (Jo et al., 2021). Water samples were transported to a laboratory installed in our accommodation using a sledge. Although the systematic effects of differences in the apparatus used for water collection on eDNA detection are not clarified, we assumed the results were comparable.

Sea-ice samples were collected from St. L on March 8, 2021. An ice block measuring 0.3 × 0.3 m in area and 0.5 m thickness was collected, and the bottom 5 cm was cut using an ice saw. This layer was composed of columnar ice, which was characterized by vertically aligned ice crystals. This layer likely formed within 2 months prior to sampling. The ice samples were placed in a clean plastic bag, transported in a cooler box to the laboratory at the Faculty of Fisheries Sciences, Hokkaido University (Hakodate, Hokkaido), and subsequently melted at 4 °C in the dark.

Seawater samples and meltwater from sea ice were filtered through a 0.45 µm pore size cartridge filter (Sterivex HV, Millipore, Billerica, MA, USA) using a peristaltic tubing pump (Masterflex 07528-10, Cole-Parmer, IL, USA) at a flow rate of 100 mL/min. Filtration stopped when the filter was severely clogged. The water volume filtered through one filter ranged from 1.3 to 8.5 L for seawater and approximately 0.3 L for meltwater from sea ice (Table 1). The filters were stored at -20 °C until DNA extraction. Prior to sample filtration, 500 mL purified water was filtered in the same manner as the negative control. All equipment used for water collection was decontaminated using a bleach solution (approximately 0.12% NaOCl) and rinsed with distilled water before use. Researchers wore new gloves and masks at each step of eDNA sampling. The workspace was wiped with a breach disinfectant spray (approximately 0.1% NaOCl) and distilled water before use.

DNA extraction, PCR, and sequencing

DNA was extracted from each filter using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the protocol of Wong et al. (2020). Detailed information on DNA extraction is provided by Kawakami et al. (2023b). DNA extract (150 µL) was obtained from each filter and stored at -20 °C.

The sequencing data used in this study were obtained from two distinct sequencing runs (RIR-005 and HFS-013, Table S1) following the Illumina paired-end index sequencing protocol on an iSeq 100 (Illumina, San Diego, CA, USA). The composition of the PCR mixture and tagging protocol for Illumina sequencing differed slightly between the runs because the methods were conducted in two distinct laboratories and the protocol had to be modified to decrease index hopping, which is the misassignment of DNA reads between samples.

As the first step, the partial mitochondrial 12S rRNA gene (approximately 170 bp) was amplified using MiFish primers (Miya et al., 2015). Selecting optimal primers for eDNA metabarcoding with adequate species coverage and taxonomic resolution is crucial for a comprehensive overview of the biodiversity at the study site (Zhu & Iwasaki, 2023). The MiFish database (v3.83, available at https://github.com/billzt/MiFish), which is a subset of GenBank sequences, completely covers the species that inhabit the Saroma-ko Lagoon, with high taxonomic coverage of 86.5% at the species level and 94.9% at the family level for fish reported in waters around Hokkaido (Amaoka et al., 2020). Therefore, using MiFish primers are expected to provide a comprehensive view of fish diversity in the Saroma-ko Lagoon.

The first PCR for RIR-005 was performed in a 20 µL reaction containing 1× PCR Buffer for KOD FX Neo (Toyobo, Osaka, Japan), 0.4 mM dNTP mix, 0.4 U KOD FX Neo polymerase, 0.29 μM of each MiFish-U primer, 0.15 μM of each MiFish-E primer, and 2 μL of the template. The PCR for HFS-013 was performed in a 12 µL reaction containing 1× PCR Buffer for KOD FX Neo (Toyobo), 0.4 mM dNTP mix, 0.4 U KOD FX Neo polymerase, 0.3 µM of each MiFish-U primers, 0.15 µM of each MiFish-E primer, and 2 µL of template. The same PCR conditions were applied: 94 °C for 2 min; 40 cycles of 98 °C for 10 s, 65 °C for 30 s, and 68 °C for 30 s, followed by a final extension at 68 °C for 5 min. Eight PCR replicates were performed for each sample. Each PCR replicate was pooled and subsequently purified using ExoSAP-IT (Thermo Fisher Scientific, Waltham, MA, USA). A no-template control (NTC) was introduced in each run of the first PCR to evaluate cross-contamination and sequencing with 2.0 μL of nuclease-free water (UltraPure; Invitrogen, Waltham, MA, USA) instead of extracted DNA.

A second PCR was performed to attach the Illumina sequence adapter and dual-index tags (i7 and i5) to the first PCR products. For RIR-005, combinational indexing was applied, with each sample having a unique combination of indices. The 13 µL reaction included 1× PCR Buffer for KOD FX Neo, 0.37 mM dNTP mix, 0.24 U KOD FX Neo polymerase, 0.27 μM of each index primer, and 2 μL of the first PCR product. For HFS-013, non-redundant indexing was applied, ensuring that each sample had completely unique indices. Although the Illumina sequencer equipped with patterned flow cell, such as iSeq 100, may generate index hopping at a high frequency (0.29% in HiSeq 2500), non-redundant indexing effectively reduces the risk (MacConaill et al., 2018). The reaction volume was 12 µL, and it contained 1× PCR Buffer for KOD FX Neo, 0.4 mM dNTP mix, 0.24 U KOD FX Neo polymerase, 0.3 µM of each index primer, and 1 µL of the first PCR product.

For each iSeq run, the second PCR products were pooled in one tube, using 2–6 µL per sample based on band brightness from gel electrophoresis, and 4 µL for the NTC and negative controls. The band brightness was visually assessed by comparing it to that of DNA ladder markers (OneSTEP Ladder 50, Nippon Gene, Tokyo, Japan). The PCR products were not diluted to normalize concentration, as DNA from high-biomass species can mask the detection of low-biomass species (Deagle et al., 2018). The pooled library was purified using AMPure XP (Beckman Coulter, Brea, CA, USA). The second PCR products with a length of approximately 370 bp, were excised using a 2% E-Gel SizeSelect Agarose Gel (Thermo Fisher Scientific) and quantified using a Qubit dsDNA HS Assay Kit and Qubit Fluorometer (Thermo Fisher Scientific). The purity was evaluated using an Agilent 2100 Bioanalyzer system (Agilent, Santa Clara, CA, USA). The library was diluted to 40 pM with Buffer EB (Qiagen) and sequenced using an iSeq 100 with a 150 bp paired-end reading and 30% PhiX spike-in (Illumina) for quality control. Additionally, negative controls for field sampling and PCR were sequenced (Table S1).

The samples were processed in a laboratory dedicated to eDNA analysis for DNA extraction and subsequent PCR preparation. PCR amplification and library preparation were conducted in a separate room to prevent cross-contamination. Researchers wore new gloves and masks at each step of molecular work. The workspace and equipment (micropipettes, incubator, and centrifuge) were wiped with DNA-OFF (Takara, Shiga, Japan) and distilled water before use. In addition, the extraction equipment (manifold, valves, and tubes) were decontaminated using a bleach solution (approximately 0.12% NaOCl), as described in the eDNA sampling section.

Sequencing data processing

Raw sequences were demultiplexed using the Generate FastQ analysis module (v.2.0.1) on iSeq and exported as Fastq files. The Fastq files were processed using the Mifish pipeline v3.91 (Zhu et al., 2023; available at: https://mitofish.aori.u-tokyo.ac.jp/mifish/) with default parameters and a 97% similarity threshold. This pipeline generates zero-radius operational taxonomic units (ZOTUs), which are valid sequences corrected for sequencing errors (Edgar, 2016) as a consequence of consecutive processes, including quality checks, assembly of sequence reads, length filtering, removal of primer sequences, denoising, and chimera removal (Zhu et al., 2023). Subsequently, this pipeline automatically assigned ZOTUs to a known species with the highest similarity based on a homology search against the implemented DNA database. ZOTUs that shared the same species names were compulsorily selected and pooled to produce a sample-by-species table. Accordingly, this assignment procedure can yield spurious results when multiple species cannot be discriminated based on sequence similarity. Therefore, to obtain robust and conservative assignments, we disregarded the results exported from the pipeline and refined the taxonomic assignments for each ZOTU as described in the taxonomic assignment section below.

To remove the probable contaminant reads, the maximum read count of each ZOTU detected in the negative controls (field blank and NTCs) was subtracted from the read count detected in the sample. Subsequently, ZOTUs were considered positive when the read count was ≥10. Prior to decontamination, read counts were corrected in proportion to the volume of the second PCR products of each sample used for library construction (Table S1). This adjustment was performed for each iSeq run, as the relationship between DNA concentration and read count could vary across runs.

Taxonomic assignment

All data handling and statistical analyses were performed using R (R Core Team, 2021) v.4.1.2, in combination with R studio (2023.12.0+369).

The ZOTUs were searched against GenBank’s nr database using BLASTN 2.14.1+ (as of 5 August 5, 2023, https://blast.ncbi.nlm.nih.gov/Blast.cgi) after combining them into a FASTA file. To reduce the risk of missing species with similar sequences owing to the dominance of sequences from a few species, the maximum number of aligned sequences displayed was set at 250. BLAST results were downloaded and parsed using a custom R script. To avoid uncertain species-associations only sequences identical to ZOTUs (100% identity) with 100% coverage were retained as candidates. Taxonomic information for these candidate species was retrieved from the National Center for Biotechnology Information (NCBI) taxonomy database using “classification” function in taxize v.0.9.100 (Chamberlain & Szocs, 2013).

To improve the accuracy of the assignment, candidate sequences derived from species found in the waters surrounding Hokkaido were selected according to a species list previously documented by Amaoka et al. (2020). ZOTUs with hit sequences from mammals were excluded from analysis. For consistency, synonyms were consolidated according to the scientific names used in the NCBI taxonomy. The ZOTUs were subsequently assigned to the lowest common ancestor (LCA), which is the lowest taxonomic level shared by the candidate species using the “lowest_common” function in the taxize. For ZOTUs annotated above the genus level, the species name was designated as “sp.” because it was unclear whether the DNA was derived from a single or multiple species. ZOTUs assigned to the same LCA collapsed to form taxonomic units unless their candidate species were completely distinct from each other.

Statistical analyses

To assess sample completeness, species accumulation curves were generated for each sample using the “iNEXT” function in iNEXT v.2.0.20 (Hsieh et al., 2016) based on corrected read counts. Given that all curves reached a plateau (Figure S1), it was assumed that the sequencing depth was sufficient to detect all the taxa contained in a filter replicate. As fluctuations in filtration volume per filter can affect the fish diversity detected by eDNA metabarcoding (Kawakami et al., 2023b), read counts were rarefied to the minimum read count within each unique sample prior to combining the results from filtration duplicates for the sample. Rarefaction was performed using the functions in vegan 2.6–6.1 (Oksanen et al., 2024): species occurrence probabilities in a rarefied community were estimated with the “drarefy” function, and species with probabilities greater than 0.5 were retained. The rarefied read count for each species was expressed as the mean of 100 iterations of the “rrarefy” function.

The taxonomic composition was compared among the samples within the under-ice seawater at St. L, between the under-ice seawater at St. L and under-ice river water at St. R, and between eDNA sources (under-ice seawater at St. L and meltwater from sea ice at St. L, under-ice river water at St. R and meltwater from sea ice at St. L). For each comparison, samples were rarefied to the minimum read count within the combination, following the same method described above. Compositional dissimilarity was evaluated using the Jaccard index (ß_jac) and Venn diagrams. ß_jac was further partitioned into turnover (β_jtu) and nestedness components (β_jne). The turnover component indicates the replacement of some species with others between samples, whereas the nestedness component indicates the extent to which the taxonomic composition of a sample with fewer species is a subset of a richer sample (Baselga, 2010). The contribution of each component to the total dissimilarity is expressed as the ratio of ß_jtu or ß_jne to ß_jac. The dissimilarity indices were calculated using the “beta.pair” function in betapart v.1.6 (Baselga et al., 2012) after converting the taxon-by-sample table into presence/absence data. Additionally, the Chao2 estimator (S_chao2), which estimates true taxonomic richness (Hsieh et al., 2016), was calculated based on results from the four samples of under-ice seawater at St. L using the “iNEXT” function (Hsieh et al., 2016) after being rarefied as previously described.

To discuss the possibility of fish overwintering within the frozen lagoon and river, the congruence of data between taxa detected by eDNA and known fish fauna was evaluated based on the National Survey on the Natural Environment, which was periodically conducted by the Ministry of the Environment, Japan, from 1973 to 2004, and the National Census on River Environments, which was continually conducted by the Ministry of Land, Infrastructure, Transport, and Tourism, Japan, after 1990. Fish species lists were downloaded from the Biodiversity Center of Japan, Ministry of the Environment (Ministry of the Environment, Japan, 1993; available at https://www.biodic.go.jp/reports2/4th/kosho/4_kosho_allm.pdf) and River Environmental Database (Ministry of Land, Infrastructure, Transport, and Tourism, Japan, 2024; available at https://www.nilim.go.jp/lab/fbg/ksnkankyo/). As no comprehensive report on fish fauna is available for the Saromabetsu River, we selected the Tokoro River as an alternative under the assumption that they have similar fish communities. In these reports, fish fauna were surveyed through fish collection, compilation of fishery catch records, visual inspection, and interviews with residents. Only taxa identified or assigned at the species level were used for comparison. Since the Actinopteri species known from these references were completely covered by the MiFish database (v.3.83), we assumed that the risk of false-negative results due to the lack of DNA data was negligible. This possibility was further evaluated by integrating environmental preferences. Habitat information of each species was retrieved from Fishbase using the “ecology” function in rfishbase v.4.1.2 (Boettiger et al., 2012) with Amaoka et al. (2020) as a complementary. The ranges of preferred parameters and absolute minimum or maximum parameter thresholds of certain species regarding water temperature and ice concentration were obtained from AquaMaps (https://www.aquamaps.org/) using aquamapsdata v.0.1.6 (Kaschner et al., 2019). As AquaMaps compiles natural occurrence data for marine species, our comparison also focused on marine fish.