Generally, the mitochondrial 16S rRNA and COI genes reflects the species diversity of animals, and COI gene has a faster base substitution rate than 16S (e.g. Palumbi 1996; Deagle et al. 2014). In brittle star, however, there is often little difference in the results of molecular phylogenetic analyses using these two regions (Boissin et al. 2008; Hunter and Halanych 2008; Cho and Shank 2010; Okanishi et al. 2018). Furthermore, since the COI gene is a coding region, many mutations accumulate in the third codon, making it relatively difficult to create universal primers to detect many species. Therefore, in this study, the mitochondrial 16S rRNA gene was selected as the target sequence for eDNA metabarcoding. This gene contains both hyper-variable and conserved regions (Gray et al. 1984; Okanishi and Fujita 2013). In order to select suitable regions of 2 primer sets for species detection by metabarcoding analysis, we referred to 1,839 16S rDNA sequences from 132 brittle star species obtained from the International Nucleotide Sequences Database (INSD). The number of species and length of the sequences were 32 with 118–160 bp for 16SOph1, and were 119 with 160–200 bp for 16SOph2 (Table 1). For developing 16SOph2, we newly sequenced partial sequences of 16S rRNA of 13 brittle star species by the same method with constructing custom reference data base (Table 1; see also “Determination of target sequences of specimens collected from the study area” below). These sequences were from species across 6 major orders and 18 superfamilies of brittle stars (Table 1). The amplified products were designed to be less than 200 bp in length, taking into account the fact that DNA in environmental water is short and fragmented (Riaz et al. 2011; Miya et al. 2015; Valentini et al. 2016). The downloaded and newly obtained sequences were aligned using Clustal W (Thompson et al. 1994) implemented in MEGA7 (Kumar et al. 2016) with a default set of parameters. In this setting, sites with gap were completely deleted from analysis. Finally, the new primer sets were designed to amplify the adjacent conserved region (16SOph1) and hyper-variable region (16SOph2) of 16S rRNA, using MEGA7 for 16SOph1 and Primer3web (Untergasser et al. 2012) for 16SOph2 with visual observation. 16SOph1 is located 5′ of 16SOph2, and 16SOph1-R partially overlaps with 16SOph2-F (Table 2; Suppl. material 4: fig. S1).

In order to test the detection power of the 2 newly designed primer sets (16SOph1, 16SOph2), a total of 1L of rearing water was obtained and mixed from individual aquarium tanks in MMBS in which 15 species of brittle stars were reared. The water was sampled from a total of 9 tanks. Eight of these tanks have a volume of ca. 2.5 L and only one larger tank has a volume of 170 L. 1L of water was collected from each of the smaller tanks (100 mL each) and from the larger tank (200 mL); as a replicate, 50 ml of water was collected twice for each 100 ml collection and 100 ml of water was collected twice for the 200 ml collection by using 300 ml measuring cups. All water samples were taken from the surface area of the tanks. Of the 8 small tanks, 1 tank contained 2 species and 2 tanks contained 3 species of brittle stars. Ophiomaza cacaotica, Astrocladus coniferus and Ophiothela danae were attached to Oxycomanthus japonicus (sea feather), unidentified sponges and unidentified octocorals, respectively. The population of each brittle star was 1 to 10 individuals or greater. As for Ophiothela danae, there were definitely more than 10 individuals, but the number was too large to determine the exact number of individuals from external observation (Suppl. material 5: table S1).

In the marine field, environmental water samples from the shallow waters of Sagami Bay facing MMBS were collected from as diverse bottom substrata as possible. In this study, we collected 1L of seawater on 7 April 2021 from each of: 1) the water surface of Moroiso (35.1558°N, 139.6050°E); and 2) its rocky bottom, ca. 2 m in depth; and 3) the muddy bottom under the pier of the MMBS (35.1576°N, 139.6121°E, ca. 2 m depth). Additionally, we also collected 5L of seawater from each of: 4) the sandy and muddy bottom depths of 84.5 m (35.1552°N, 139.5766°E) on 24 September 2021; and 5) 250 m (35.1164°N, 139.5706°E) and 270 m (35.1164°N, 139.5701°E) on 8 July 2021, respectively, using the research vessel (R/V Rinkai-Maru) of MMBS with a Niskin water sampling bottle of 5L volume (K Engineering Co. Ltd.) (Figs 1, 2). The reason we sampled 5L at these depths was because we expected that deep water would be cleaner and contain less DNA than shallow water (e.g. McClenaghan et al. 2020). These are areas where the fauna of brittle stars has been well studied (e.g. Matsumoto 1917; Irimura 1982; Fujita et al. 2006).

Figure 1.

A–C. Map of Japan (A), showing the location of the 2 sea water sampling sites of Misaki Marina Biological Station (B) and 2 sea water sampling sites at shallow waters and aquarium tank (C) (Koajiro, Misaki, Miura, Kanagawa Prefecture).

Figure 2.

A–D. Natural sea water sampling operations at the Pier with a bucket, ca. 3 m (A); at water surface (B) and under water (C) of Moroiso with a 1 L polychlorinated bottle, ca. 2 m; at water bottom of Jogashima, ca. 80 m, 250 m and 270 m with 5L water sampling bottle (D). E, F. Aquarium tanks at Misaki Marine Biological Station, whole view (E) and an individual tank (F).

Environmental water at the pier was collected using 12 L buckets, and when sampling under the pier with the buckets, a 10 m rope was fastened to a bucket and 500 mL of bottom water was collected twice to gather 1L of seawater (Fig. 2A). Sampling at the rocky bottom at a depth of 2 m and at the water surface at Moroiso was performed by scuba diving, and environmental water was collected underwater using 1L polychlorinated bottles (Fig. 2B, C). This “1L” volume follows a previous study (Komai et al. 2019) that did similar research on the Japanese coast. These buckets, polychlorinated bottles, and measuring cups were decontaminated thoroughly with sodium hypochlorite solution (final concentration approximately 1%) before use. The 5L Niskin water sampler bottle was decontaminated with fresh water before use (Fig. 2D). After sampling, each water sample was stirred well in the field with benzalkonium chloride solution (final concentration 0.01%) to suppress the degradation of eDNA (Yamanaka et al. 2017) and sent refrigerated to the laboratory at Environmental Research & Solution Co. Ltd. on the same day, for subsequent analysis. For the water sampling at 80 m, 250 m and 270 m, 5L of Milli-Q water was brought on the boat from the laboratory in a clean bottle as a blank sample, and the Milli-Q water was transferred to the sampling bag at the field. The blank samples were treated in the same way as field samples in subsequent analyses.

The water sample was filtered through a glass fiber filter, Whatman GF/F filters (GE Healthcare, Japan, average pore size 0.7 um). DNA extraction and purification were performed using the DNeasy and Tissue kit (Qiagen) following the protocol of Tsuji et al. (2020). After extraction, DNA was stored at -25 °C. eDNA extracted from all the environmental seawaters was subjected to the first-round PCR (1^st PCR) and the second-round PCR (2^nd PCR) in order to append amplified sequences with 3 kinds of adaptor sequences: 1) primer-binding sites for sequencing, 2) dual-index sequences to distinguish amplicons, and 3) sequences for binding to the flowcells of the Illumina MiSeq (Illumina, CA, USA). Before performing the 1^st PCR, the DNA extract was repurified using magnetic beads using the protocol of the manufacturer of Ampure XP. The 1^st PCR was carried out with a 12 μl reaction volume containing 6.0 μl KAPA HiFi HotStart ReadyMix (KAPA Biosystems, MA, USA), 1.4 μl of each 16SOph primer (5 μM primer F/R), 1.2 μl sterile Milli-Q water and 2.0 μl eDNA template. The thermal cycle profile after an initial 3 min denaturation at 95 °C was as follows: 40 cycles of denaturation at 98 °C for 20 s, annealing at 60 °C (for 16SOph1)/56 °C (for 16SOph2) for 15 s and extension at 72 °C for 15 s with the final extension at the same temperature for 5 min. The 1^st PCR products were purified using Magnetic Beads (AMPure XP, Beckman) in order to remove dimers and monomers following the manufacturer’s protocol. Subsequently, the purified products were quantified using QuantiFluor ONE dsDNA System (Promega), and diluted to 0.1 ng/μl using Milli-Q water, and the diluted products were used as a template for the 2^nd PCR. The 2^nd PCR was conducted with a 12 μl reaction volume containing 6.0 μl KAPA HiFi HotStart ReadyMix, 0.83 μl each primer (5 μM), 2.34 μl sterile distilled H₂O and 2.0 μl template. The thermal cycle profile after an initial 3 min denaturation at 95 °C was as follows: 12 cycles of denaturation at 98 °C for 20 s, annealing and extension combined at 72 °C for 15 s, with the final extension at the same temperature for 5 min. In order to monitor contamination during the process of PCRs, blank samples were prepared. During the 1^st PCR, a filtration blank (FB), an extraction blank (EB) and a PCR blank (1B) with 2.0 μl Milli-Q water instead of template eDNA were added. All the libraries containing the target region and the 3 adapter sequences were mixed in equal volume and the pooled libraries were size-selected from approximately 230 bp for 16SOph1 and 320 bp for 16SOph2 (including variable region+primer annealing site+NNNNNN+2^nd primer binding site+slow cell binding site) using a 2% E-Gel Size Select agarose gel (Invitrogen, CA, USA). The distribution of the products of the size-selected libraries was measured and confirmed using a BioAnalyzer (Agilent Technologies) and the concentration of the size-selected libraries was measured using a QuantFlour ONE dsDNA System, and they were sequenced on the MiSeq platform using a MiSeq Reagent v2 Kit for 2×150bp PE or 2×250 PE (Illumina, CA, USA) following the manufacturer’s protocol.

Equipment for filtration was immersed in 5% sodium hypochlorite solution, washed with tap water, and further washed with Milli-Q water. Gamma-sterilized filter microtips were used for DNA extraction, and the DNA purification process was left to an automated extraction system (QIA cube, QIAGEN) to eliminate human error as much as possible. 1^st PCR to amplify the DNA of the target taxa and the post-PCR processing were performed in separate rooms. To avoid as much as possible “Tag jumps”, in which the wrong combination of tags is used in sequencing by next generation sequencing (Schnell et al. 2015), forward and reverse indexed primers (DualIndex) were used to identify samples pooled after 2^nd PCR. Each pair of F primer and R primer was designed to differ by at least two bases.

All data preprocessing and analysis of MiSeq raw reads were performed using USEARCH v10.0.240 (Edgar 2010) according to the following steps. 1) Both forward and reverse reads were merged by aligning them using the fastq_mergepairs command. During this process, low-quality tail reads with a cut-off threshold set at a quality (Phred) score of 2, too short reads (< 50 bp) after tail trimming and those paired reads with too many differences (> 5 positions) in the aligned region (ca. 130 bp for sequences from aquarium tank by 16SOph1 and ca. 160 bp and ca. 210 bp for other sequences by 16SOph1 and 16SOph2, respectively) were discarded; 2) primer sequences were removed from those merged reads using the fastx_truncate command; 3) those reads without the primer sequences underwent quality filtering using the fastq_filter command to remove low quality reads with an expected error rate of > 1% and too short reads of < 50 bp; 4) the preprocessed reads were dereplicated using the fastx_uniques command and all singletons, doubletons and tripletons were removed from the subsequent analysis following the recommendation of Edgar (2010); 5) the dereplicated reads were denoised using the unoise3 command and all putatively chimeric and erroneous sequences were separated from the subsequent amplicon sequence variant (ASV) assignment; 6) finally, all processed reads with a sequence identity of > 80% with our “custom database 2+3” were assigned to sequences of brittle stars (given the registered name of the database) using the usearch_global command (Fig. 3; see below for the explanation of taxonomic assignment).

Figure 3.

Graphical summary of the dataset construction protocol used in this study.

As reference sequences for taxonomic assignment, we used a custom database including a total of 8,975,113 mitochondrial 16S rRNA gene sequences (https://doi.org/10.6084/m9.figshare.22139909.v1; “custom database 1”) from Metazoa, Plant, Fungi, Eubacteria and Archaea and 1,340 sequences from 33 known brittle star families (Suppl. material 1; “custom database 2”), which were downloaded from INSD on 17 December 2021. Therefore, custom database 1 includes the sequences of custom database 2. In addition to those published sequences, we independently determined 16SOph1 sequences from 56 species and 16SOph2 sequences from 44 species (Table 4) and added those sequences to the reference database (68 reference sequences in total; morphologically identified to 59 species from 33 genera of 16 families of 4 orders; “custom database 3”). These sequences were newly obtained by Sanger sequencing in this study (see also “Determination of target sequences of specimens collected from the study area” below for detailed methods).

A phylogenetic tree was then reconstructed using the “sequences of brittle stars” from each study area to confirm its monophyly. All sequence sets of brittle stars obtained in the data processing from each water sample were aligned using the Clustal W algorithm in MEGA7. All missing bases were scored as gaps and the sites with gaps were completely excluded from phylogenetic analysis. The substitution models were computed for each sequence set with the “find best-fit model of nucleotide substitution” option by MEGA7. Maximum likelihood analysis (ML) with 1000 bootstrap replicates was performed with MEGA7 to reconstruct the ML phylogenetic tree. The trees were visualized with MEGA 7. Monophyly of each node in the phylogenetic trees was considered to be supported if bootstrap was higher than 80%. Node bootstrap values lower than 79% were considered as not monophyletic (Suppl. material 4: figs S2–S10).

To determine the species “threshold distance” we measured the genetic distance for each 16SOph1 and 16SOph2 region of custom database 3 to obtain the minimum genetic distance interspecies. We refrain from using custom database 2 since it was expected to include far fewer sequences from Sagami Bay, unidentified sequences such as “Ophiuroidea sp.”, and misidentified sequences. If a genetic distance between a sequence from custom database 3 and a sequence from the eDNA was estimated to be within the “threshold distance”, the latter sequence was given the taxon name in the database, and treated as “species with certain identification”. Taxa including all biota other than the brittle stars were searched with the 2 primer sets as sequences that matched our custom database 1 with more than 97% similarity.

The DNA sequences of the 16S rRNA gene were determined for 59 morphologically identified species that had been collected in Sagami Bay and from similar environments in the seas around Japan during the last 10 years. This included specimens from areas where water samples were collected (including tanks) for the present study. The method of DNA extraction followed that of Okanishi and Fujita (2013). The primer set 16Sar (5′-CGC CTG TTT ATC AAA AAC AT-3′) and 16Sbr (5′-CCG GTC TGA ACT CAG ATC ACG T-3′) (approximately 500 bp on the 5′ side of the 16S rRNA gene) was used and its amplicon region covered the target regions of both 16SOph1 and 16SOph2 (Palumbi 1996). For specimens that 16Sar and 16Sbr would not amplify, an alternative newly designed internal primer (Oph-16S-R1) covering 16SOph1 and 16SOph2 was used. The primer sets used for amplifying the sequence including 16SOph1 and/or 16SOph2 are summarized in Table 3 and Suppl. material 4: fig. S1. The optimum cycling parameters for those primer sets consisted of an initial denaturation step of 95 °C/2 min followed by 41 cycles of 95 °C/30 s, 48 °C/40 s and 72 °C/1 m with final extension step at 72 °C/10 m, which was followed by storage at 10 °C (Okanishi and Fujita 2013). The PCR products were separated from excess primers and oligonucleotides using Exo-SAP-IT (GE Healthcare), following the manufacturer’s protocol. Sequences amplified with 16Sar and 16Sbr were sequenced bidirectionally and sequence products were run on a 3730xI DNA Analyzer (Thermo Fisher Scientific). The accession numbers of DNA Data Bank of Japan (DDBJ) of the examined specimens were deposited in INSD with collection locality data (Table 4). For all examined materials, coordinates, depth, date, collector, gear, and identifier were provided in as much detail as possible (Suppl. material 3).

As a result of aligning the mitochondrial 16S rRNA sequences from datasets of 132 species of brittle stars (Table 1), we were able to design new primer sets for the conserved region of about 110 bp and the hyper-variable region of about 180 bp. These primer sets were named “16SOph1” and “16SOph2”, respectively (Table 3). Melting temperatures and G/C contents of 16SOph1-F/16SOph1-R/16SOph2-F/16SOph2-R were 63.8 °C/62.5 °C/59.9 °C/64.3 °C and 52%/58%/61%/53%, respectively. The length of the 16SOph sequence data obtained from environmental water was 111–115 bp (112 bp on average) for 16SOph1 and 180–195 bp (185 bp on average) for 16SOph2. These sequences were registered in DDBJ/EMBL/NCBI (Table 4).

The genetic distance of custom database 2 was examined, and 22,360 (16SOph1) and 6,510 (16SOph2) sequence combinations were calculated with distance of 0, which indicates they matched 100% in sequence. Among these, many combinations were found that differed at the order level, indicating that they were probably misidentified (Suppl. material 5: tables S2, S3). The genetic distance of the sequences in custom database 3 showed that most combinations with distance 0 were the same species. However, the genetic distance between Amphiura koreae and Amphioplus rhadinobrachius was also calculated as 0 (16SOph1; Suppl. material 5: table S4). In addition, even within the same species (Ophiothela danae), genetic distances of 1.87% (16SOph1) and 0.63% (16SOph2) were calculated depending on the difference of collection sites, and some degree of genetic distance was calculated between A. trachydisca and A. koreae (1.87% for 16SOph1; 1.2% for 16SOph2), and between A. rhadinobrachius and A. trachydisca (1.87% for 16SOph1) (Suppl. material 5: tables S4, S5). Therefore, in this study, these genetic distances (1.87% for 16SOph1 and 1.2% for 16SOph2) were considered to indicate intraspecific variation, and the three species of Amphiuridae (A. koreae, A. trachydisca and A. rhadinobrachius) were considered synonyms (see also Discussion for this taxonomic issue).

Except for these three species of Amphiuridae, the smallest genetic distances were 2.86% (16SOph1) and 6.15% (16SOph2) between Stegophiura sladeni and Stegophiura sterea (Suppl. material 5: tables S4, S5). Therefore, in this study, we set 97.14% and 93.85% as the threshold values for the species boundary of 16SOph1 and 16SOph2, respectively, and the eDNA sequences in the database with similarities exceeding this value were considered the same species.

After merged, quality-filtered, dereplicated and denoised the raw reads, we obtained, respectively for 16SOph1/16SOph2, a total of 362,300/574,409 usable processed reads for aquarium tank, 141,629/152,491 for the sea surface of Moroiso, 142,856/157,812 for the sea bottom of Moroiso, 134,055/281,242 for the Pier, and 360,480/160,237 for the 80 m depth, respectively. The number of usable processed reads obtained for 250–270 m sample was 177,844 (only 16SOph1).

We compared the processed reads (hereafter “reads”) obtained from eDNA of environmental waters (including aquarium water tanks) at MMBS with our custom database 2+3. As a result, we recorded the following number of sequences of brittle stars from this database that matched with >80% similarity (hereafter “sequences”) the eDNA reads: a total of 20 sequences from the tank (comprising a total of 338,367 reads, the range was 167–278,769, and the average was 21,147 reads), 37 sequences from Moroiso (totals of 94,056 and 109,680 reads, the range was 162–25,618 and 62–26,844, and the average was 2,542 and 2,964 reads from surface and bottom, respectively), 19 sequences from the Pier (a total of 86,056 reads, the range was 15–30,422, and the average was 7171 reads), 52 sequences from 80 m depth (a total of 182,951 reads, the range was 4–108,554, and the average was 4,691 reads), and 20 sequences from 250 m and 270 m depths (a total of 37,211 reads, the range was 5–10,824, and the average was 1,162 reads) (Suppl. material 5: tables S6–S10). The sampling points at 250 m and 270 m depths were geographically close to each other, and we combined the data from the 2 areas.

We constructed maximum likelihood phylogenetic trees for these sequences of brittle stars in each study area and filtered the obtained monophyletic sequences by the “threshold” value, resulting in: 16 phylogenetic clades in the tank, 14 of which could be considered as species with certain identification (Suppl. material 4: fig. S2; Suppl. material 5: table S6); 18 phylogenetic clades in Moroiso, 14 of which could be considered as species with certain identification (Suppl. material 4: fig. S3; Suppl. material 5: table S7); 12 phylogenetic clades in Pier, 9 of which could be considered as species with certain identification (Suppl. material 4: fig. S4; Suppl. material 5: table S8); 39 phylogenetic clades in the 80 m sample, 18 of which could be considered as species with certain identification (Suppl. material 4: fig. S5; Suppl. material 5: table S9); 16 phylogenetic clades in the 250 and 270 m sample, 10 of which could be considered as species with certain identification (Suppl. material 4: fig. S6; Suppl. material 5: table S10). In Moroiso, for the clade of Macrophiothrix longipeda, the bootstrap value was 67. However, considering that the similarity range of the clade was 97.3–100, we considered this clade to be a species with certain identification. In the Pier sample, for the clade of Amphiura iridoides the bootstrap value was 79. However, because the difference in nucleotide sequence was within 2 bp, we consider it to be a species with certain identification. In the 80 m depth sample, the range of identity of 3 sequences of Amphiura digitula obtained was 93.8–100%, and this value is below the 16SOph1 threshold. However, these 3 sequences form a monophyletic clade in the phylogenetic tree, and their maximum similarity was 100%. Therefore, they were considered to be a species with certain identification (Suppl. material 5: tables S6–S10).

Among them, 10/13 species, 8/15 species, 3/9 species, 7/18 species, and 5/10 species were confirmed to be reared in the tank, or distributed in Moroiso, the Pier, 80 m, and 250–270 m, respectively (Suppl. material 5: tables S6–S10).

We compared the processed reads obtained from eDNA of environmental waters (including the aquarium tanks at MMBS) with our custom database 2+3. As a result, a total of 36 sequences of brittle stars from the tank (comprising a total of 522,242 reads, the range was 18–482,876, and the average was 34,816 reads), 42 sequences from Moroiso (totals of 6,147 and 13,460 reads, the ranges were 13–2,397 and 8–4,353, and the averages were 409 and 897 reads from the surface and bottom, respectively), 14 sequences from the Pier (a total of 11,186 reads, the range was 8–4,065, and the average was 1,016 reads), and 16 sequences from 80 m (a total of 755 reads, the range was 4–575, and the average was 75 reads) (Suppl. material 5: tables S11–S13) were found to match our custom database with more than 80% similarity. Our maximum likelihood phylogenetic tree based on sequences in each study area and filtration of the obtained monophyletic sequences by the “threshold” value resulted in: 15 phylogenetic clades in the tank, 15 of which could be considered as species with certain identification (Suppl. material 4: fig. S7); 15 phylogenetic clades in Moroiso, 15 of which could be considered as species with certain identification (Suppl. material 4: fig. S8); 11 phylogenetic clades in the Pier, 11 of which could be considered as species with certain identification (Suppl. material 4: fig. S9); and 10 phylogenetic clades in the 80 m depth sample, 10 of which could be considered as species with certain identification (Suppl. material 4: fig. S10).

Among them, 11/15 species, 10/15 species, 4/11 species, and 1/10 species were confirmed to be reared in the tank, or distributed in Moroiso, at the Pier, and at 80 m depth, respectively (Suppl. material 5: tables S11–S14). In 80 m, the range of identity of 7 sequences of Amphichilus trichoides was 92.1–100%, and this value is below the 16sOph2 threshold. However, these 7 sequences form a monophyletic clade in the phylogenetic tree, and their maximum similarity was 100%. Therefore, they were treated as a certainly identified species.

The results of the search of all other biota other than the brittle stars for 16SOph1 indicated that almost all of the detected species were brittle stars in all environmental waters, and 1 insect species was detected at Moroiso and in the deep sea, respectively, which may have been due to contamination during the experimental analysis (Fig. 4). Two species of nemerteans were detected in the aquarium tanks. This may have been due to the species that were reared together with the brittle stars in the tank. On the other hand, only brittle stars were detected at the Pier and 80 m (Fig. 4). Non-metazoa biota (Plant, Fungi, Eubacteria, and Archaea) were not detected using 16SOph1.

Figure 4.

Compositions of OTUs (species) and reads of the all biota (including detected brittle star species with certain identification in this study) detected from eDNA around MMBS by 16SOph1 and 16SOph2.

The proportion of OTUs (species)/reads putatively assigned to 1) brittle star was 87.5%/99.9% (Tanks), 91.6%/99.9% (Pier), 100%/100% (Moroiso surface), 90.9%/99.9% (Moroiso bottom), 100%/100% (80 m) and 91.6%/99.9% (deep sea); 2) Echinodermata other than brittle star was 0% for all sites; 3) other metazoans was 12.5%/>0.001% (Tanks), 8.3%/>0.001% (Moroiso surface), 9%/>0.001% (Moroiso bottom) and 8.3%/>0.001% (deep sea) (Fig. 4; Table 5). A wide variety of species were detected with 16SOph2. In addition to the brittle stars: 20 species (8 phyla, 11 classes) were detected in the tank, 43 species (6 phyla, 9 classes) in the surface of Moroiso, 43 species (8 phyla, 11 classes) in the bottom of Moroiso, and 46 species (6 phyla, 10 classes) in the Pier. On the other hand, only 1 species of crustacean (Arthropoda) was detected at 80 m (Fig. 4). The number of brittle star species detected was more than one-quarter of that in the tank, but less than one-quarter of that in Moroiso and at the Pier, while the number of brittle star species detected was more than 90% of that at 80 m (Fig. 4; Table 5). Non-metazoan biota was detected in the tank (Proteobacteria), Moroiso surface (Pyramimonadophyceae and Mamiellophyceae), Moroiso bottom (Pyramimonadophyceae), and Pier (Pyramimonadophyceae).

The proportion of OTUs (species)/reads putatively assigned to 1) brittle star was 40.5%/99.4% (Tanks), 22.4%/40.6% (Pier), 18.9%/20.1% (Moroiso surface), 22.5%/16.1% (Moroiso bottom), 90.9%/100% (80 m), and 91.6%/98.3% (deep sea), respectively; 2) those assigned to Echinodermata other than brittle star were 10.8%/0.003% (Tanks), 13.7%/2.1% (Pier), 6.8%/1.3% (Moroiso surface), 6.4%/0.001% (Moroiso bottom), and 0%/0% (80 m), respectively; 3) other metazoans were 43.2%/0.001% (Tanks), 60.3%/56.8% (Pier), 72.4%/78.3% (Moroiso surface), 64.6%/83.5% (Moroiso bottom), and 9%/1.6% (80 m), respectively; 4) non-metazoans were 5.4%/>0.001% (Tank), 6.4%/>0.001% (Moroiso surface), 3.4%/0.001% (Moroiso bottom), and 8.3%/>0.001% (80 m) (Fig. 4; Table 5).

A total of 38 taxa of brittle stars with precise species-level identification (without “Ophiuroidea sp.: RCMBA668.2231” which is registered as an unidentified species in INSD) were obtained from water eDNA from the shallow (including aquarium tank water eDNA) to the deep sea waters around the Misaki Marine Biological Station. The species detected with the primer pairs 16SOph1 and 16SOph2 were respectively: 14/15 species detected in the aquarium water, 10/14 species in the water surface of Moroiso, 11/13 species on the sea bottom of Moroiso, 9/11 species on the pier, and 18/10 species at a depth of 80 m (Table 4; Suppl. material 5: tables S6–S14; Suppl. material 4: figs S11–S15). These results indicate that the 16SOph2 primer has higher detection power for species in shallow water than the 16SOph1 primer. On the other hand, at a depth of 80 m, the number of species detected by 16SOph1 was more than double the number detected by 16SOph2. Based on these results, eDNA metabarcoding was performed using 16SOph1 at 250–270 m depth, which is a sandy and muddy substrate similar to that at 80 m depth. As a result, we obtained 10 brittle star species (including one species of “Ophiuroidea sp.”) from these depths that were accurately identified (Table 4; Suppl. material 4: fig. S16).

The difference between the 2 primer sets was clearly shown in the taxa detected other than the brittle stars. In all of the water eDNA samples analyzed in this study, more than 90% of the species and reads detected by 16SOph1 were brittle stars. Considering that the sequences of insects obtained from Moroiso and the deep sea were contaminants, the only other organisms’ sequences that were surely picked up were those of 2 nemertean species that were reared together with the brittle stars in the tanks (Fig. 4; Table 5).

On the other hand, among the species detected with 16SOph2, the proportion of brittle stars in shallow water was less than half in all environments, and less than one-quarter in the Moroiso and the Pier. This is the same as the trend in obtained reads: in Pier and Moroiso, the number of reads of brittle stars was less than half. In Tank, however, most of the reads were brittle stars. On the other hand, by detection with16SOph2, mollusca, bryozoans, and crustaceans rather than other echinoderms (asteroids, crinoids, echinoids, and holothuroids) were more abundant in Moroiso and Pier. The proportions of non-metazoans (plants and bacteria) in both 16SOph1 and 16SOph2 were quite small in all environments (Fig. 4; Table 5).

These results suggest that 16SOph1 performs better at detecting mainly brittle stars, while 16SOph2 is able to detect a broader range of marine metazoa. In addition to the detected brittle stars, other species were found to be distributed in the Misaki area, suggesting that this is not a result of contamination. However, since this study focuses on brittle stars, we do not discuss the results of other taxa in detail here.

In addition, at a depth of 80 m, even with 16SOph2, only 1 species of crustacean was detected. One hypothesis to explain this is that there is a very large amount of eDNA from brittle stars at the bottom of the deeper sea area, much more than eDNA from crustaceans. In order to verify this hypothesis, we need to collect water samples in many more deep-sea environments and test other primers in the future.

In the last 10 years, 59 morphologically identified species have been collected from the MMBS (Table 4). As a result of the eDNA metabarcoding analysis in this study, we were able to detect at least 38 of these species, which means that we have succeeded in developing a method for obtaining biodiversity information without collecting specimens or conducting visual observations. There have been few studies on eDNA metabarcoding of deep-sea organisms at depths deeper than 200 m based on comparison with detailed collected specimen data. Therefore, this study provides a method to expand our knowledge of environmental monitoring in the ocean, which contains more than 90% of the Earth’s biosphere.

It should also be noted that the genetic distance analysis of sequences obtained by using 16SOph1 and 16SOph2 suggested that Amphioplus rhadinobrachius, Amphiura koreae, and Amphiura trachydisca should be synonyms. The taxonomy of Amphiuridae, to which the three species belong, has been confusing, and in fact, the recent phylogenomic analysis of brittle stars showed that the genus-level phylogeny including Amphiura and Amphioplus is polyphyletic (O’Hara et al. 2017). The results of our study are therefore not surprising. We here refrain from synonymising these three species since it is beyond the scope of our study. The species-genus level classification of Amphiuridae should be taxonomically reviewed in the future.

The sequences with 80% < identity < 97% with our custom database 2+3 were obtained for 16SOph1 as follows: 2 species at the tank, 4 species at Moroiso, 3 species at the pier, 20 species at a depth of 80 m, and 6 species at depths of 250–270 m. Because of the ambiguity of identification, we refrain from mentioning these species in this study. However, considering the number of the detected species, it is suggested that there are many brittle stars that have not yet been discovered from Sagami Bay. On the other hand, the sequences with identity less than 93.85% (uncertain identification) were not detected by 16SOph2.

In terms of the total number of reads (and their percentages compared to all detected metazoan taxa): 338,367 (99.9%)/522,016 (99.4%) in the aquarium tank; 109,680 (99.9%)/5,791 (16.1%) in Moroiso (including surface and sea bottom), 86,056 (100%)/11,186 (20.1%) in the Pier, and 182,951 (100%)/755 (98.3%), were detected using 16SOph1 and 16Soph2, respectively (Table 5). Considering these results, it can be seen that although the number of detected reads of Oph16S2 was high in the aquarium tank, the overall trend was for those of Oph16S1 to be higher. The proportions of reads of brittle stars among all metazoans in the tank and at 80 m were almost the same for both primers, but those in Moroiso and the Pier were much higher for 16SOph1 than for 16SOph2 (Table 5). These results clearly showed the characteristics of 16SOph1 targeting the conserved region and 16SOph2 targeting the hyper-variable region. The number of species with certain identification was higher for 16SOph2, but the number of species with uncertain identification was higher for 16SOph1 (Suppl. material 5: tables S6–S14). This may be attributed to the shorter target sequence of the former primer, which could more easily pick up DNA fragmented in environmental water. This is also consistent with the fact that 16SOph1 had a higher number of detected reads. However, the fact that 16SOph2 amplified the DNA sequences of other marine animals other than brittle stars was unexpected. This may be a result of chance. In addition, when phylogenetic trees were constructed, 16SOph2 tended to have longer branches between each monophyletic clade than 16SOph1 and higher bootstrap values supporting the monophyly of each clade (Suppl. material 4: figs S2–S10). Thus, 16SOph1 and 16SOph2 have their own advantages and disadvantages. Therefore, it is expected that more accurate measurement of the marine environment will be possible by using one primer or the other for each environment, or by using both. In the following sections, we would like to discuss the detected species as follows. For the taxonomic notes for the detected species, see Suppl. material 2.

In this study, we examined the detection power of the developed primers from the aquarium waters where the brittle stars were actually reared (Fig. 2E, F). As a result, 12 (16SOph1) and 13 (16SOph2) taxa were detected by each primer among 15 reared species (Suppl. material 4: fig. S11). Ten species were common to both primers (Suppl. material 5: tables S6, S11). Ophiactis macrolepidota1 was detected only by 16SOph1, and Stegophiura sterea, Ophiactis profundi, and Ophioplocus japonicus were detected only by 16SOph2. Of these, O. japonicus, for example, is a common species that occurs in the shallow waters around MMBS. However, considering the fact that this species was detected only by 16SOph2 at Moroiso and the Pier, it is clear that there was a difference in the detection trend for each primer.

Ten (16SOp1) and 14 species (16SOp2) were detected in the surface water of Moroiso (Suppl. material 5: tables S8, S9). Other than Amphiura sinicola, the species detected by 16SOph1 were also detected by 16SOph2, and Ophiuroglypha kinbergi, Ophiothela danae, Ophiocentrus tokiokai, and Ophiactis lymani were detected only by 16SOph2. Therefore, 16SOph2 was more powerful in detecting the brittle stars in the water surface of Moroiso (Suppl. material 4: figs S12, S13; Suppl. material 5: tables S9, S10).

Eleven species (16SOph1) and 13 species (16SOph2) were detected on the sea bottom of Moroiso (Suppl. material 5: tables S8, S9). Among them, Ophioplocus japonicus, Ophiocentrus tokiokai, Ophiactis lymani, and Ophiothrix exigua were detected only by 16SOph2, and Ophiomaza cacaotica and Ophionereis dubia were detected only by 16SOph1. Also on the sea bottom of Moroiso, 16SOph2 was more powerful in detection of taxa of brittle stars than 16SOph1 (Suppl. material 4: fig. S13; Suppl. material 5: tables S8, S9).

Nine taxa (16SOph1) and 11 taxa (16SOph2) were detected from the Pier (Suppl. material 5: tables S10, S11). Ophiactis lymani, Amphichilus trichoides, and Macrophiothrix longipeda were detected only by 16Soph2, and Ophiomastix mixta was detected only by 16SOph1. The detection power of 16SOph2 was higher than that of 16SOph1 (Table 4; Suppl. material 4: fig. S14; Suppl. material 5: tables S10, S11).

At a depth of 80 m, 18 species (16SOph1) and 10 species (16SOph2) of brittle stars were detected. In this area, species detected by 16SOph2 were generally also detected by 16SOph1, except for Ophiothrix nereidina, which was only detected by 16SOph2. Therefore, the detection power of 16SOph1 was higher than that of 16SOph2 in this area (Suppl. material 4: fig. S15; Suppl. material 5: tables S12, S13).

Nine species were detected from 250–270 m depth (Suppl. material 5: table S14). This sampling area is also an area frequently surveyed by dredges of the laboratory research boat, and at least 23 species (considering synonyms) have been recorded there so far. In this study, Ophiothrix nereidina and Ophiosparte gigas were newly detected by eDNA metabarcoding (Suppl. material 4: fig. S16; Suppl. material 5: table S14).

False negatives (not detected even though present) and false positives (detected even though absent) are important issues for eDNA metabarcoding because they can lead to underestimation or overestimation of biodiversity (e.g. Komai et al. 2019). The results of the present analysis showed that many of the species inhabiting the shallow waters were actually detected (Suppl. material 5: tables S6–S14). On the other hand, many of the brittle star species recorded from 80 m and 250–270 m were not detected in the present study; and conversely, many sequences of species that are supposed not to exist or not to be collected in these area were detected (Table 4; Suppl. material 5: tables S9, S13, S14; Suppl. material 2).

Factors causing false negatives include 1) DNA of the target species was not collected during sampling, and 2) DNA of the target species was not amplified well by PCR (e.g. Komai et al. 2019). Regarding the first factor, the difficulty of collecting water using a water sampler can be mentioned. Because the water sampler is lowered to the seafloor on a single string, the deeper the water depth, the greater the amplitude of the water sampler near the seafloor. Therefore, in deep water, the slight rocking of the boat is greatly transmitted to the water sampler on the seafloor, and consequently the sampler tends to move away from the seafloor. In this study, we did not observe the actual water sampling, so we do not know if we were able to drop the water sampler just above the seafloor. In the future, we expect to be able to collect more eDNA by collecting a large amount of water as close to the seafloor as possible using an ROV and larger water sampler.

In order to solve the second problem, it is necessary to develop optimal experimental conditions using 16SOph primers. In the future, increasing the number of PCR replicates and using multiple annealing temperatures, as well as the revision of primer sequences, will allow us to detect more target eDNA (Komai et al. 2019).

Finally, a third factor, DNA shedding, must be taken into account. In seawater, DNA from all organisms is not always present. An organism’s DNA is shed at some level, and the rate of shedding should differ among taxa. Therefore, no matter how efficient the water sampling is or how suitable the primers are, some organisms that live in special environments and have specific ecologies will remain undetectable. In the future, it will still be important to supplement data by collecting organisms that inhabit the study area, rather than relying solely on environmental DNA for monitoring.

One of the most important points for accurate assessment of biodiversity by eDNA metabarcoding is to improve the reference sequence database, which is essential for taxonomic assignment. In this study, we created a custom database for the preparation of primers by adding our own sequence data to the INSD database. Currently, 350 species of brittle stars are known from Japanese waters (Okanishi 2016; Okanishi and Fujita 2018a, b, c; Okanishi et al. 2019; Okanishi et al. 2020a, b; Okanishi and Kohtsuka 2021), which is about 17% of the currently known brittle star species (32 families, ca. 2100 species; Stöhr et al. 2012). In addition, there are many undescribed cryptic species or species complexes, such as Amphipholis squamata (Boissin et al. 2008) and Asteronyx loveni (Okanishi et al. 2018), and Ophiactis savignyi (Roy and Sponer 2002), and such species were also detected in this study (e.g. Ophiothrix panchyenchyta1 and 2; Table 4; Suppl. material 2). In order to solve the cryptic species problem, a thorough quantitative morphometrical approach and re-examination of type specimens with observation of internal ossicle morphology are required, along with DNA phylogenetic analysis as in a study by dos Santos Alitto et al. (2019).

Overrepresentation of a few frequently sequenced species in INSD is also an issue that should be noted. In this study, more than 95% of the sequences that were identified with certainty were compared to sequences that we obtained in this study or from a previous study (e.g. Okanishi and Fujita 2013). Therefore, we do not anticipate that our study was affected by the overrepresented sequences accumulated on the INSD due to primer bias. However, in future studies, it will be necessary to take into account the effects of such sequences and remove the overrepresented sequences appropriately.

In this study, we sequenced the partial 16S rRNA genes of 59 brittle star species (three of these species are suspected to be synonyms; Suppl. material 2) obtained from around the study area, among which 18 species (31.5%) did not match eDNA (Table 4). Brittle stars in Sagami Bay have been studied for more than 100 years, starting with Matsumoto’s research (e.g. Matsumoto 1917; Irimura 1982; Fujita et al. 2006). However, the fact that many unknown sequences (comprising at least 33 species) were obtained by eDNA metabarcoding indicates that the taxonomic study of brittle stars is still incomplete even in this region.

There are many taxa, such as the brittle star, that have not received much attention so far, but are, in fact, rich in DNA and have great potential as research targets for environmental monitoring. Until now, much data on COI sequences has been accumulated in animals as so-called “DNA barcoding regions” (Hebert et al. 2003). In recent years, however, it has been recognized that mitochondrial rRNA sequences are suitable for eDNA monitoring in a variety of organisms (e.g., Miya et al. 2015; Komai et al. 2019). For example, in fish, information on the mitochondrial genome was originally accumulated for each species. Therefore, it was possible to convert the data into reference sequence data (Iwasaki et al. 2013; Miya et al. 2015). For other taxa, however, information on the entire mitochondrial genome has not yet been accumulated, and more time will be needed to construct a reference sequence database for eDNA monitoring.

In addition to the conventional COI-based DNA barcoding project, it will be necessary to determine the data of mitochondrial genes, including 16S rRNA, for each species to optimize the usefulness of eDNA metabarcoding analysis. This will require a significant taxonomic update by taxonomists.

In this study, we constructed a database based on specimens actually collected and organized by taxonomists, and laid the foundation for an eDNA metabarcoding method for brittle star in the deep sea. The oceans occupy most of the Earth’s biosphere, and marine animals account for most of the Earth’s animal biomass (Bar-On et al. 2018). Therefore, in order to accurately measure the biodiversity of the Earth, it is essential to accurately understand the diversity of marine organisms. This study is expected to be a pioneering study in deep-sea research, where eDNA metabarcoding has rarely been performed due to the difficulty of access, despite the fact that the volume of the deep sea is the largest in the ocean.