All seawater sampling at the six sites was conducted in compliance with Japanese laws and regulations as well as local ones.

As part of the various outreach activities at the Natural History Museum and Institute, Chiba, the present study recruited six pairs of parents and children on the website (http://www2.chiba-muse.or.jp/NATURAL/) to collect eDNA samples from marine environments in Chiba and the surrounding prefectures. A seminar and a workshop were conducted on July 4, 2021, to provide an overview on eDNA studies and to secure high-quality eDNA samples. A sampling kit was distributed to each group during the workshop and its use was demonstrated. Moreover, examples of appropriate or inappropriate sites for water sampling in coastal areas were demonstrated, but no specific regions or locations were mentioned to the participants. Thus, they freely chose water sampling sites, depending on their personal preferences of residence and summer holiday locales. All relevant collection data are summarized in Suppl. material 1: Table S1.

For seawater collection, low-tech bucket sampling was utilized, comprising various small commercial buckets fastened to a rope (provided by the participants). Prior to seawater sampling, the participants wore disposable gloves on both hands and assembled a set of on-site filtration kits consisting of a Sterivex filter cartridge (pore size 0.45 µm; Merck Millipore, MA, USA) and a 50-ml disposable syringe with a Luer lock connector (TERUMO, Tokyo, Japan). Subsequently, they fixed the end of the rope fastened to the bucket and collected surface seawater by casting and retrieving the bucket filled with seawater.

The participants performed on-site filtration using the above kit (filter cartridge + syringe) to collect and concentrate the eDNA on the filter membrane inside the cartridge. They removed the filter cartridge from the syringe and drew approximately 50-ml of seawater into the syringe by pulling the plunger, reattaching the filter cartridge to the syringe, and pushing the plunger for filtration of the seawater. This step was repeated until the final filtration volume reached 1000 mL. When the filter was clogged before reaching 1000 mL filtration, the total volume of filtered seawater was recorded.

After on-site filtration, an outlet port of the filter cartridge was closed with a 3-mm diameter rubber cap (KOKUGO, Tokyo, Japan), filled with 1.6-mL RNAlater (Thermo Fisher Scientific, DE, USA) into the cartridge from an inlet port of the cartridge using a disposable capillary pipette (AS ONE, Tokyo, Japan) to prevent eDNA degradation, and an inlet port was closed with a screw cap (TERUMO, Tokyo, Japan) for preservation. They transported the filtered cartridges to their homes in a portable cooler with ice packs and kept them at 4 °C in their fridges prior to shipment. Subsequently, they returned the cartridges to the museum using a refrigerated courier, and the six cartridges were stored at –20 °C until eDNA extraction was performed.

For more details of the eDNA collection, see “Environmental DNA sampling and experiment manual Version 2.1 (The eDNA Society 2019)

Details of the laboratory protocol can be found in the Suppl. materials 5: Supplementary methods. The workspace and equipment were thoroughly sterilized prior to all laboratory experiments. Filtered pipette tips were used and eDNA-extraction, pre- and post-PCR manipulations were conducted in three different dedicated rooms that were physically separated from each other to safeguard against cross-contamination from PCR products.

eDNA was extracted from the filter cartridges using a DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) following the methods developed by Miya et al. (2016) with slight modifications. After removing the preservative in the cartridge, the filter cartridge was lysed using proteinase K. The collected DNA extract (~900 μL) was purified using the DNeasy Blood and Tissue kit following the manufacturer’s protocol and the final elution volume was set to 200 µL. An extraction blank (EB) was also created during this process.

This study employed two-step PCR for paired-end library preparation using the MiSeq platform (Illumina, CA, USA). In general, we followed the methods developed by Miya et al. (2015) and subsequently modified by Miya and Sado (2019). For the first round of PCR (1^st PCR), a mixture of the following six primers was used: MiFish-U-forward, reverse, MiFish-E-forward-v2, reverse-v2, MiFish-U2-forward, and reverse. These primer pairs amplify a hypervariable region of the mitochondrial 12S rRNA gene (ca. 172 bp; hereafter called “MiFish sequence”) and append primer-binding sites (5´ ends of the sequences before six Ns) for sequencing at both ends of the amplicon.

The 1^st PCR was performed with eight technical replicates for the same eDNA template to minimize PCR dropouts. A 1^st PCR blank (1B) was also prepared during this process, in addition to EB. After completing the 1^st PCR, an equal volume of PCR products was pooled from each of the eight replicates, and the pooled products were purified, quantified, and diluted to 0.1 ng/µL using Milli Q water, and the diluted products were used as templates for the second round of PCR (2^nd PCR). For the two blanks (EB, 1B), the 1^st PCR products were purified in the same manner.

The 2^nd PCR was performed to append dual-index sequences and flow cell binding sites for the MiSeq platform. A 2^nd PCR blank (2B) was also prepared during this process, in addition to EB and 1B. The six PCR products as well as the three blank samples from the 2^nd PCR products were pooled along with other samples from different projects. Subsequently, the pooled dual-indexed libraries were electrophoresed on an agarose gel, and the target amplicons (~370 bp) were excised. The concentration of the size-selected libraries was measured, diluted to 10.0 pM, and sequenced on the MiSeq platform using a MiSeq v2 Reagent Kit for 2 × 150 bp PE (Illumina, CA, USA) following the manufacturer’s protocol.

All raw DNA sequence data and associated information were deposited in the DDBJ/EMBL/GenBank database and are available under accession number DRA012840.

Data preprocessing and analysis of raw MiSeq reads from the MiSeq run were performed using PMiFish ver. 2.4 (https://github.com/rogotoh/PMiFish.git; Miya et al. 2020). Forward (R1) and reverse (R2) reads were merged while discarding low-quality tail reads with a cut-off threshold set at a quality (Phred) score of 2, too short reads (<100 bp) after tail trimming, and those paired reads with too many differences (>5 positions) in the aligned region (ca. 65 bp). Primer sequences were removed from those merged reads and those reads without the primer sequences underwent quality filtering to remove low quality reads. The preprocessed reads were dereplicated, and all singletons, doubletons and tripletons were removed from the subsequent analyses to avoid false positives (Edgar 2010). The dereplicated reads were denoised to generate amplicon sequence variants (ASVs) that removed all putatively chimeric and erroneous sequences (Callahan et al. 2017).

The ASVs were subjected to taxon assignments to species names (molecular operational taxonomic units; MOTUs) with a sequence identity of >98.5% with the reference sequences (two nucleotide differences allowed) and a query coverage of ≥90%. An incomplete reference database necessitates this clustering step, which enables the detection of multiple MOTUs for identical species names. Such multiple MOTUs were annotated with “gotu1, 2, 3…” and all of these outputs (MOTUs plus U98.5 MOTUs) were tabulated with read abundances. ASVs with sequence identities of <80% (saved as “no hit”) were excluded from the above taxon assignments and downstream analyses because all of them were found to be non-fish organisms. MiFish DB ver. 43 was used for taxon assignment, comprising 7973 species distributed across 464 families and 2675 genera.

To refine the above taxon assignments, family level phylogenies were reproduced from MiFish sequences from MOTUs and reference sequences (contained in the MiFish DB ver. 43) belonging to these families. For each family, representative sequences (most abundant reads) from MOTUs were assembled, all reference sequences were added from that family, and saved in FASTA format. The combined FASTA-formatted sequences were subjected to multiple alignments using MAFFT 7 (Katoh and Standley 2013) with a default set of parameters. A neighbor-joining (NJ) tree was subsequently constructed with the aligned sequences in MEGA X (Stecher et al. 2020) using Kimura two-parameter distances.

A total of 67 family-level trees were visually inspected and taxon assignments were revised in the following manner. For those U98.5 MOTUs placed within a monophyletic group consisting of a single genus, the unidentified MOTUs were named after that genus, followed by “sp.” with sequential numbers (e.g., Pagrus sp. 1, sp. 2, sp. 3...). For the remaining MOTUs ambiguously placed in the family-level tree, the unidentified MOTUs were named after that family, followed by “sp.” with sequential numbers (e.g., Sparidae sp. 1, sp. 2, sp. 3...).

The final list of detected species can be found in Suppl. material 2: Table S2.

The differences in fish community structures were visualized using non-metric multidimensional scaling (NMDS) with 999 separate runs of real data. For NMDS, community dissimilarity was calculated based on incidence-based Jaccard indices, and NMDS stress was used to confirm the representation of NMDS ordination. The “metaMDS” function of “vegan” ver. 2.5.6 package (Oksanen et al. 2019) was used for NMDS.

To evaluate whether the detected fish communities reflected the marine environments of the six sampling sites (Fig. 1), three types of salinity preferences (salt, salt/brackish, salt/brackish/freshwater) were assigned for each detected species using information available from FishBase (Froese and Pauly 2019). A PMiFish pipeline (ver. 2.4) automatically retrieved the relevant information (water area, salinity preference, depth range) through FishBase in a default setting using “rfishbase” (Boettiger et al. 2012).

Figure 1.

A) the location of the survey area and schematic flow paths of the Kuroshio (red) and Oyashio (blue) currents; B) the location of the six survey sites (Sts. 1–6) in Boso and Miura peninsulas. Map data 2021 Apple Inc.

To evaluate the biogeographic characteristics of the detected fish communities from the six sampling sites, the center of the geographic distribution for each detected species was calculated by averaging the latitudes of the northern and southern limits in the Northern Hemisphere following Masuda (2008), with reference to the distribution data in Nakabo (2013). Those data can be found in Suppl. material 3: Table S3.

As a baseline dataset, a faunal inventory was compiled for coastal marine fish from the region including the survey areas based on multiple museum collections and literature surveys (Suppl. material 4: Table S4). Five collections were selected from neighboring natural history museums to represent coastal marine fish fauna in Chiba Prefecture: Natural History Museum and Institute, Chiba (CBM) and its Coastal Branch (CNMH), National Museum of Nature and Science, Tokyo (NSMT), Kanagawa Prefectural Museum of Natural History (KPM), and Yokosuka City Museum (YCM). Subsequently, those fishes collected from the coastal marine areas of Chiba prefecture were extracted from databases of these collections using the search word “Chiba,” wherein fishes inhabiting freshwaters only were manually excluded (Suppl. material 4: Table S4). For the literature surveys, eight core references were designated (Miya et al. 1994a, b, 1995; Hagiwara and Kimura 2006; Kohno et al. 2011; Nakabo 2013; Aoki et al. 2016; Shitamitsu et al. 2019).

Six groups of parents and children were selected from 22 groups of applicants, such that their places of residence were scattered. Those six groups freely chose their sampling sites (Sts. 1–6 in Fig. 1), performing seawater sampling and on-site filtration during the summer holiday of 2021, from July 26 to August 21 (Suppl. material 1: Table S1).

St. 1 is located at the southern end of Kujukuri Beach – a shallow, vast sandy beach with a total length of 66 km along the Pacific coast – and is strongly affected by the intrusion of the cold Oyashio water along the coastline (Yang et al. 1993). Although the seabed of St. 1 is entirely sandy, wave-dissipating blocks protect the shore from coastal erosion (Fig. 2). Moreover, the seawater is often turbid due to the inflow of muddy water from an estuary of the Isumi River approximately 1 km south (Fig. 2). The filtered cartridge appeared darkest (Fig. 2) and clogged at a filtration volume of 600 mL

Figure 2.

Enlarged maps of the six survey sites (Sts. 1–6), their snapshots, and corresponding filtered cartridges. The participants took photos from Sts. 1, 3, 5, and 6 on the survey date (Suppl. material 2: Table S2), and those of Sts. 2 and 4 were taken by one of us (MM) on September 28, 2021. Map data 2021 Apple Inc.

Sts. 2–5 are located in the southern tips of the Boso and Miura peninsulas, and are greatly affected by the warm Kuroshio current flowing northward along the Pacific coast as well as its branches (Soh 2003; Fig. 1). In contrast to St. 1, these sites are located within or surrounded by rocky shores with intricate coastlines (Fig. 2), and the seawaters are generally clear, even within fishing ports (Sts. 2 and 4). The filtered cartridges remained relatively clear and were not clogged, with a maximum filtration volume of 1000 mL (Fig. 2).

St. 6 is located at the mouth of the Hanami River, deep within Tokyo Bay (Figs 1, 2). The innermost part of Tokyo Bay was originally a vast tidal flat, and the natural coast has disappeared because of repeated reclamation and revetment for over 350 years (Endoh 2004), transforming into a large industrial area in the last 50 years (Fig. 2). The filtered cartridge appeared darker, comparable to St. 1 (Fig. 2), and became clogged at a filtration volume of 600 mL.

Based on the eDNA extracts, MiFish fragments were successfully amplified in the 1^st PCR, with distinct, putatively fish bands being observed around ca. 310 bp (Miya et al. 2020) for all the six samples. Dual-index sequences and flowcell-binding sites were appended in the 2^nd PCR and the target amplicons (ca. 370 bp) were purified and quantified for the six libraries. MiSeq paired-end sequencing (2 × 150 bp) of the six libraries – in conjunction with an additional 29 libraries (total = 35) – yielded a total of 6,135,029 reads, with an average of 97.2% base calls with Phred quality scores of ≥30.0 (Q30; error rate = 0.1% or base call accuracy = 99.9%). This run was highly successful, considering that the manufacturer’s guidelines (Illumina Publication no. 770-2011-001 as of May 27, 2014) are >80% bases ≥Q30 at 2 × 150 bp.

Of the 6,135,092 reads, 531,306 were assigned to the six libraries, and the number of raw reads for each library ranged from 67,618 to 123,743, with an average of 88,551 reads. After merging two overlapping paired-end fastq files [654,497 reads (98.5%)], the primer-trimmed sequences were subjected to quality filtering to remove low-quality reads [650,461 reads (97.9%)]. The remaining reads were dereplicated for subsequent analysis and single to tripletons were removed from the unique sequences. Thereafter, the reads were denoised to remove putatively erroneous and chimeric sequences. The remaining 567,980 reads (85.5% of the raw reads) were subjected to taxon assignments after rarefaction to the minimum number of reads (67,681). Of these, 531,306 reads (93.5% of the denoised reads) were putatively considered fish sequences, while BLAST searches indicated that non-fish sequences [36,674 reads (6.5%)] primarily consisted of mammals (i.e., cows, pigs, and humans) and a few unknown sequences. The three negative controls (EB, 1 B, and 2 B) were subjected to the same analysis pipeline, whereby they did not yield any denoised reads.

Following the automatic taxon assignments, the family level NJ trees were visually inspected and the species names in the list were revised. The final list included 258 detections, assigned to 140 species across 66 families and 118 genera (Suppl. material 2: Table S2). All detected species can be found in the fish inventory comprising 948 species across 158 families and 493 genera (Suppl. material 4: Table S4), with the exception of unnamed species (family or generic names + sp.). The latter may be attributed to the lack of reference sequences or congeneric species with little or no interspecific sequence variations in the MiFish fragments, such as the rockfish genus Sebastes (Sebastidae), some flying fishes (Exocoetidae), the surfperch genus Ditremma (Embiotocidae), and the pufferfish genus Takifugu (Tetraodontidae) (Miya et al. 2020).

The number of species detected at the six sampling sites ranged from 14 (St. 1) to 66 (St. 5), with a mean of 43 (Suppl. material 2: Table S2). The number at the four sites along the rocky shores (49–66 spp. at Sts. 2–5) exceeded the two sites at a sandy shore (14 spp. at St. 1) and an inner bay (19 spp. at St. 6) (Fig. 3). For convenience, the rocky shores, sandy shore, and inner bay are hereafter called “area(s)” individually or collectively.

Figure 3.

Scatter plot of the numbers of fish species detected from the fish eDNA metabarcoding. The numbers were not corrected for the filtered volumes (1000 mL for Sts. 2–5 and 600 mL for Sts. 1, 2) and are plotted against latitudes. Orange, blue, and green dots represent the data from sandy-shore, rocky-shore, and inner-bay areas, respectively. Note that the rocky-shore areas are mostly located along the southern coastline, while sandy-shore and inner-bay areas are located in the northern coastlines of Boso Peninsula. Thus, the gradient of the species richness does not necessarily relate to latitudes.

Of the 140 detected species, only Blackhead Seabream (Acanthopagrus schlegelii) and Grass Puffer (Takifugu niphobles; putatively including few congeneric species owing to low interspecific genetic sequence variations; Miya et al. 2020) were observed across the three areas. The former species is one of the most popular targets for anglers in Japan, while the latter is a small, ubiquitous poisonous fish that is not of much importance to anglers (https://fishing-fishing.com/2020/01/22/5783/). Approximately half of the detected species at the sandy-shore and inner-bay sites (7 of 14 spp. at St. 1 and 10 of 19 spp. at St. 6) were also identified in other areas, whereas 89% of the detected species at the rocky shore sites (110 of 123 spp. at St. 2–6) were unique to this area (Fig. 4). Within the latter area (rocky shore), 13–17 species (22–33%) were unique at each site, and the rest of the species were also found at two to four sites in this area. Furthermore, NMDS ordination analysis indicated that dissimilarities in the fish communities of these three areas were well separated in two-dimensional space (Fig. 5; NMDS stress = 0).

Figure 4.

Venn diagram of the number of detected species from sandy-shore, rocky-shore, and inner-bay areas. Total numbers of the detected species for each area are indicated within parentheses.

Figure 5.

Non-metric multi-dimensional scaling (NMDS) for fish communities detected in four sites along rocky-shore (Sts. 2–5, blue), sandy-shore (St. 1, green) and inner-bay (St. 6, orange) areas. NMDS stress ≈ 0.

As expected from the geographical locations of the sampling sites (Fig. 2), fish tolerant of freshwater and saltwater accounted for approximately one-third of Sts.1 and 6 (Fig. 6). Of these fish species populating multiple habitats, Flathead Mullet (Mugil cephalus cephalus) and Yellowfin Goby (Acanthogobius flavimanus) appeared in common, although the marine environments of the two stations were significantly different (eroded sandy shore facing open ocean near the estuary vs. river mouth in the highly urbanized bay area). Other marine and brackish water species detected at the two stations (Suppl. material 2: Table S2) were frequently found in museum collections and references (Suppl. materials 2, 4: Tables S2, S4). In contrast to these two sites, marine fish were predominant at Sts. 2–5, occupying >90% of the component species (Fig. 3). Most of them were commonly found in rocky shores along the southern coast of Boso and Miura peninsulas (Miya et al. 1994a, b, Hagiwara and Kimura 2006; Shitamitsu et al. 2019) (Suppl. material 2, 4: Tables S2, S4).

Figure 6.

Proportions of the three categories of salinity preferences for detected fish species from each study site. Salinity-preference data were taken from FishBase (Froese and Pauly 2019) using “rfishbase” (Boettiger et al. 2012).

The center of the geographic distribution of each component species was calculated from the average latitude of the northern and southern limits from literature records (Suppl. material 3: Table S3). The distribution patterns have been summarized in a box plot for each station (Fig. 7). The average latitudes of the component fish species from St. 1 – which are strongly affected by the cold Oyashio water – are biased to the north, while those from Sts. 2–5, which are strongly affected by the warm Kuroshio Current, were slightly biased to the south. Such a trend is ambiguous in the highly urbanized site (St.6), despite its northernmost location.

Figure 7.

Box plots depicting variations in the center of geographic distribution for each detected species, calculated from an average of the latitudes of the southern and northern limits in the Northern hemisphere provided in Nakabo (2013). Orange, blue, and green bars represent the data from sandy-shore, rocky-shore, and inner-bay areas, respectively. Total numbers of the detected species for each area are indicated within parentheses.