Quantitative assessment of multiple fish species around artificial reefs combining environmental DNA metabarcoding and acoustic survey

Since the early 1970s, many artificial reefs (ARs) have been deployed in Japanese coastal waters to create fisheries grounds. Recently, researchers began to use environmental DNA (eDNA) methods for biodiversity monitoring of aquatic species. A metabarcoding approach using internal standard DNAs [i.e., quantitative MiSeq sequencing (qMiSeq)] makes it possible to monitor eDNA concentrations of multiple species simultaneously. This method can improve the efficiency of monitoring AR effects on fishes. Our study investigated distributions of marine fishes at ARs and surrounding stations in the open oceanographic environment of Tateyama Bay, central Japan, using qMiSeq and echo sounder survey. Using the qMiSeq with 12S primers, we found higher quantities of fish eDNAs at the ARs than at surrounding stations and different fish species compositions between them. Comparisons with echo sounder survey also showed positive correlations between fish eDNA concentration and echo intensity, which indicated a highly localized signal of eDNA at each sampling station. These results suggest that qMiSeq is a promising technique to complement conventional methods to monitor distributions of multiple fish species.

www.nature.com/scientificreports/ Also, this method requires less field time and is more likely to detect the presence of target species than the other survey methods 6,9 . Recently, the metabarcoding approach, using high-throughput sequencing technologies (e.g., Illumina MiSeq) and universal primer sets, has been utilized to reveal the fish diversity in a given area 6,10 . MiFish, a universal polymerase chain reaction (PCR) primer set on the mitochondrial 12S rRNA gene, is especially popular for fish eDNA metabarcoding and has detected more fish species than other competing primers 11,12 . Although the eDNA metabarcoding approach has become a cost-and labor-effective approach for estimating aquatic biodiversity, it has one limitation, i.e., the quantity of eDNA cannot be estimated from the number of eDNA sequence reads obtained by high-throughput sequencing partly due to several experimental problems such as PCR inhibitions 7 . However, some studies have proposed an approach using an internal standard DNA to overcome this limitation. Ushio et al. 7 tested whether additions of internal standard DNAs (i.e., known amounts of short DNA fragments from fish species that have never been observed in a sampling area) to eDNA samples can create sample-specific regression lines between the numbers of sequence reads and DNA copies. Based on regression lines for each sample, they converted the sequence reads of DNAs from unknown fish species (i.e., non-standard fish species) to the number of fish DNA copies in each sample. This study showed that the estimated numbers of total and two dominant fish (Japanese anchovy and jack mackerel) DNA copies were significantly positively correlated with the numbers of DNA copies quantified by qPCR. Therefore, this method [i.e., quantitative MiSeq sequencing (qMiSeq)] can contribute to the quantitative assessment of multispecies fish eDNA. Some studies have shown that the eDNA quantity could be a proxy for the abundance or biomass of macro-organisms, not only in closed systems like tank and mesocosm 13 , but also in marine ecosystems 9,14 . The qMiSeq with MiFish primer could be one of the most cost-and labor-effective approaches for the quantitative monitoring of multiple fish species in marine and coastal ecosystems, including ARs.
Here, we investigated marine fish distributions around ARs in Tateyama Bay, central Japan, using eDNA metabarcoding, catch record of fishery, and echo sounder method. The qMiSeq using internal standard DNAs and MiFish primers was applied to assess whether the distribution of multiple fish species around ARs could be monitored based on the spatial variation in the estimated quantities of eDNA (i.e., identification of fish species and quantification of the number of fish eDNA copies simultaneously). Specifically, the following were tested: (1) whether fish communities detected by eDNA metabarcoding match with the catch record of the coastal fishery, (2) whether the DNA quantities of multiple fish species calculated by the qMiSeq decrease with increasing distance from an AR, 3) whether the spatial variation of fish DNA quantities correlates to the relative distribution estimated by the acoustic survey.

Materials and methods
Study site, field survey, and in situ filtration. The field survey was performed in Tateyama Bay (34° 60′ N, 139° 48′ E), central Japan, in the proximity of the Kuroshio warm current facing the Pacific Ocean (Fig. 1). This area has many artificial reefs (ARs) created to improve fishing efficiency for fishers. Among the ARs, we focused on one high-rise steel AR (AR1), with a height of 30 m, where fish tended to aggregate ( Fig. 1 and S1). Sampling stations were set up at the AR1 and at six linear distant points extending northeast and southwest. These stations were named E150, E500, E750, W150, W500, and W750, where "W" or "E" and the number of each station name represented northeast or southwest and distance in meters from the AR1, respectively (Table S1 and Fig. 1). Another station was set up at a second AR (AR2: 25 m height) 220 m from AR1 because we found AR2 by chance during the survey (Table S1 and Fig. 1), and it might affect the eDNA concentration at other stations.
We conducted water sampling at eight stations for eDNA analysis and performed an acoustic survey for estimating relative fish density using research vessel Takamaru (Japan Fisheries Research and Education Agency: FRA) on May 23, 2018. We started the echo sounder survey at the eastern part of the bay and continued it during the water sampling ( Fig. 1). Although the echo sounder survey could not differentiate between fish species, we collected this data to assess the association between the estimated concentration of fish eDNA and the echo intensity measured by the echo sounder. Water sampling began at E750, then continued along the transect line to E150, AR1, W150, W500, W750, before going back to AR2. At each sampling station, we collected 10 L of seawater from both the middle and bottom layers by one cast of two Niskin water samplers (5L × 2 samples) and measured vertical profiles of water temperature and salinity with a conductivity-temperature-depth sensor (RINKO profiler, JFE Advantech Co., Ltd.). We subsampled 2L seawater from the 5 L seawater of Niskin sampler using measuring bottle and remaining 3 L seawater was used for pre-wash of measuring bottle and filtration devices. Two 2L samples were collected from two Niskin water samples, and then immediately filtered using a combination of Sterivex filter cartridges (nominal pore size = 0.45 μm; Merck Millipore) through an aspirator (i.e., the two filters were subsets of a single water collection) in a laboratory on the research vessel. After filtration (average time of 15 min), an outlet port of the filter cartridge was sealed with an outlet luer cap, 1.5 ml RNAlater (Thermo Fisher Scientific Inc., Waltham, MA) was injected into the cartridge using a filtered pipette tip to prevent eDNA degradation, and an inlet port was also sealed with an inlet luer cap 14 . The Niskin water samplers were bleached before each water collection using a commercial bleach solution while filtering devices (i.e., filter funnels and measuring cups used for filtration) were bleached after every filtration. We filtered 2L MilliQ water with a filter funnel and measuring cup as a field negative control to test for possible contamination. The filter cartridges were placed in a freezer immediately after filtration until eDNA extraction. In total we collected and filtered 32 eDNA samples (eight stations × two depth layers × two replicates). Disposable latex or nitrile gloves were worn during sampling and replaced between each sampling station.
DNA extraction and purification. Workspaces were sterilized prior to DNA extraction using 10% commercial bleach, and filter tip pipettes were used to safeguard against cross-contamination. Following the method  15 , the eDNA was extracted and purified. Briefly, after removing RNAlater inside the cartridge using a centrifuge, proteinase-K solution was injected into the cartridge from the inlet port, and the port was re-capped with the inlet lure cap. The eDNA captured on the filter membrane was extracted by constant www.nature.com/scientificreports/ stirring of the cartridge at a speed of 20 rpm using a roller shaker placed in an incubator heated at 56 °C for 20 min. The eDNA extracts were transferred to a 2-ml tube from the inlet of the filter cartridges by centrifugation. The collected DNA was purified using a DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol. After the purification, DNA was eluted using 100 μl of the elution buffer (buffer AE). All DNA extracts were frozen at − 20 °C until paired-end library preparation.

Preparation of internal standard DNAs. Five artificially designed and synthetic internal standard
DNAs, which were similar but not identical to the region of any existing fish mitochondrial 12S rRNA, were included in the library preparation process to estimate the number of fish DNA copies [i.e., quantitative MiSeq sequencing (qMiseq)] 7,16 . They were designed to have the MiFish primer-binding regions as those of known existing fishes and to have the conserved regions in the insert region. Variable regions in the insert region were replaced with random bases so that no known existing fish sequences had the same sequences as the standard sequences. The standard DNA size distribution of the library was estimated using an Agilent 2100 BioAnalyzer (Agilent, Santa Clara, CA, USA), and the concentration of double-stranded DNA of the library was quantified using a Qubit dsDNA HS assay kit and a Qubit fluorometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). Based on the quantification values obtained using the Qubit fluorometer, the copy number of the standard DNAs was adjusted as follows: Std. A (100 copies/µl), Std. B (50 copies/µl), Std. C (25 copies/µl), Std. D (12.5 copies/µl) and Std. E (2.5 copies/µl). Then, these standard DNAs were mixed.
Paired-end library preparation. Two PCR-level negative controls (i.e., each with and without internal . These primer pairs co-amplify a hypervariable region of the fish mitochondrial 12S rRNA gene (around 172 bp) and append primer-binding sites (5′ ends of the sequences before five Ns) for sequencing at both ends of the amplicon. The five random bases were used to enhance cluster separation on the flow cells during initial base call calibrations on the MiSeq platform. The thermal cycle profile after an initial 3 min denaturation at 95 • C was as follows (35 cycles): denaturation at 98 • C for 20 s; annealing at 65 • C for 15 s; and extension at 72 • C for 15 s, with a final extension at the same temperature for 5 min. Eight replications were performed for the 1st PCR, and the replicates were pooled to minimize the PCR dropouts. The 1st PCR products from the eight tubes were pooled in a single 1.5-ml tube. Then, we sent the 1st PCR products to IDEA consultants, Inc. to outsource the following MiSeq sequencing processes. The pooled products were purified and size-selected for 200-400 bp using a SPRIselect (Beckman Coulter, Inc.) to remove dimers and monomers following the manufacturer's protocol. The second-round PCR (2nd PCR) was carried out with a 24 µl reaction volume containing 12 µl of 2 × KAPA HiFi HotStart ReadyMix, 2.8 µl of each primer (5 µM), 4.4 µl of sterilized distilled H 2 O, and 2.0 µl of template. We used the following two primers to append the dual-index sequences (8 nucleotides indicated by Xs) and flowcell-binding sites for the MiSeq platform (5′ ends of the sequences before eight Xs): 2nd-PCR-forward (5′-AAT GAT ACG GCG ACC ACC GAG ATC TAC ACX XXX XXX XAC ACT CTT TCC CTA CAC GAC GCT CTT CCG ATC T-3′); and 2nd-PCR-reverse (5′-CAA GCA GAA GAC GGC ATA CGA GAT XXX XXX XXG TGA CTG GAG TTC AGA CGT GTG CTC TTC CGA TCT-3′). The thermal cycle profile after an initial 3 min denaturation at 95 • C was as follows (12 cycles): denaturation at 98 • C for 20 s; combined annealing and extension at 72 • C for 15 s, with a final extension at 72 • C for 5 min. The concentration of each second PCR product was measured by quantitative PCR using TB Green Fast qPCR Mix (Takara inc.). Each sample was diluted to a fixed concentration and combined (i.e., one pooled 2nd PCR product that included all samples). The pooled 2nd PCR product was size-selected to approximately 370 bp using BluePippin (Sage Science). The size-selected library was purified using the Agencourt AMPure XP beads, adjusted to 4 nM by quantitative PCR using TB Green Fast qPCR Mix (Takara Bio Inc.), and sequenced on the MiSeq platform using a MiSeq v2 Reagent Kit (2 × 150 bp) (Illumina, Inc.).
Data preprocessing and taxonomic assignment. The raw MiSeq data were converted into FASTQ files using the bcl2fastq program provided by Illumina (bcl2fastq v2.18). The FASTQ files were then demultiplexed using the command implemented in Claident 17 . We adopted this process rather than using FASTQ files demultiplexed by the Illumina MiSeq default program in order to remove sequences with low-quality scores and PCR artifacts (chimeras).
The processed reads were subjected to a BLASTN search against the full NCBI database. We excluded unique sequences of the following settings: the sequence belonged to organisms other than bony fishes, sharks, and rays; the sequence similarity between queries and the top BLASTN hit was < 98.5%; the sequence length was less than ≤ 150 bp. After BLASTN searches, assembled sequences assigned to the same species were clustered, and we considered the clustered sequences as operational taxonomic units (OTUs). www.nature.com/scientificreports/ After this automatic taxonomic assignment, we modified the taxonomic assignments because the program may assign a single species name to an OTU, whereas closely related species cannot be distinguished using the 12S rRNA gene 6 . Following Yamamoto et al. 6 , we assigned genus or a higher taxonomic rank to an OTU if the OTU was composed of sequences of different species. For example, some query sequences were aligned to two closely related species, Scomber japonicus and S. australasicus. In this case, we assigned Scomber spp. to this OTU. However, even such ambiguously assigned OTUs were re-assigned to a single species if only one species from the taxa occurs around Tateyama Bay. For example, Acanthopagrus schlegelii and A. sivicolus cannot be distinguished using the 12S rRNA gene, but we were able to assign the species rank Acanthopagrus schlegelii to this OTU because A. schlegelii is the only species from the genus Acanthopagrus in Tateyama Bay.
To remove possible contaminants, we removed the OTUs whose sequence reads were < 0.05% of the total reads as a noise for each sample. Therefore, singletons were also eliminated by this process. Three species (Oncorhynchus nerka, Gadus chalcogrammus, and Trachurus trachurus) were detected in this study but they cannot obviously exist in Tateyama Bay and are commonly consumed as food items by humans. We considered them as contaminations and excluded them as well.
Catch record of coastal fishery. To confirm the feasibility of 12S metabarcoding, we compared the fish community of eDNA metabarcoding with a catch record of a net set near the study site (Fig. 1). The detection rate was defined as the proportion of species detected by metabarcoding compared to the species collected by the net. Because our eDNA survey was conducted in May 2018, we used catch-record data from the same period (May 2018).
Estimation of DNA copy numbers. According to Ushio et al. 7 and Ushio 15 , the procedure to estimate DNA copy numbers consisted of two parts: (i) linear regression analysis to examine the relationship between sequence reads and the copy numbers of the internal standard DNAs for each sample and (ii) the conversion of sequence reads of non-standard fish DNAs to estimate the copy numbers using the result of the linear regression for each sample.
Linear regressions were used to examine how many sequence reads were generated from one fish DNA copy through the library preparation process. Note that a linear regression analysis between sequence reads and standard DNAs was performed for each sample, and the intercept was set as zero. The regression equation was: MiSeq sequence reads = regression slope × the number of standard DNA copies [/µl]. The number of linear regressions performed was thirty-three (= the number of fish DNA samples + field negative controls with standard DNAs), and thus thirty-three regression slopes were estimated in total.
The sequence reads of non-standard fish DNAs were converted to copy numbers using sample-specific regression slopes estimated by the above regression analysis. The number of non-standard fish DNA copies was estimated by dividing the number of MiSeq sequence reads by a sample-specific regression slope (i.e., the number of DNA copies = MiSeq sequence reads/regression slope). A previous study demonstrated that these procedures provide a reasonable estimation of DNA copy numbers using high-throughput sequencing 7 .
Acoustic survey of fish school using echo sounder. We estimated relative fish density using a quantitative echo sounder following the acoustic survey of Yamamoto et al. 9 . We used the echo sounder, KFC-6000 (Sonic Co., Tokyo, Japan), with a transducer (frequency, 38 kHz; beam type, split-beam; beam width, approximately 7°; pulse duration, 0.3 ms; ping rate, 0.4 s). The transducer was set at the bottom of the research vessel at a depth of 2.1 m to avoid cavitation bubbles generated by the research vessel. The acoustic devices were operated during the entire survey cruise, and all signals were recorded. The average ship speed was ~ 10 knots between sampling stations, although the ship slowed down when approaching a sampling station and completely stopped when water samples were being collected.
We eliminated noise from the obtained echo intensity data using Echoview ver. 11.0 (Echo-view Software Pty. Ltd., Tasmania, Australia). Signals between the sea bottom and 1.0 m above the sea bottom were eliminated to avoid possible confounding with the acoustic dead zone. In addition, the signals of Niskin water samplers and RINKO profiler were also eliminated. Finally, we eliminated the signals of bubbles generated by the movement of the screw propeller. After eliminating these noises, we applied a threshold of −70 dB based on visual inspection of the echograms to remove background noise 18 and assumed that the echo intensity of fish was more than this threshold. This threshold was also utilized in other studies 19,20 . Then, we extracted and averaged the echo intensity data for ten minutes after the time of each eDNA sampling and used it to assess the association with eDNA concentration. Note that the obtained echo intensity data would include echo signals from a variety of fish species. Data analyses. We performed three types of univariate analysis using generalized linear models (GLMs).
First, we used likelihood ratio tests to assess variation in the number of OTUs, eDNA copy numbers of total fish, demersal fish, pelagic fish, and five dominant OTUs, as well as echo intensity among sampling stations. When the effects of sampling stations were significant, the Tukey-HSD test was used to determine the significant differences between sampling stations. The second analysis assessed the number of OTUs and eDNA copy numbers of total, demersal, pelagic, and five dominant OTUs in relation to distance from the AR1, sampled depth layer (i.e., middle or bottom layers, represented as 'depth'), and the interaction between the distance and depth layer. Although our sampling stations contained two ARs, we focused on AR1 in the second analysis because more fish aggregated around AR1 and transport of eDNAs from AR1 may have more impact than those from AR2. The third analysis examined the eDNA copy numbers (each eDNA copy number of total fish, demersal fish, pelagic fish, or five dominant OTUs described below) in relation to echo intensity (S V ), sampled depth layer, and the interaction between the echo intensity and depth layer. Because echo intensity can correlate to distance from the We also performed multivariate analyses to determine differences in the fish community structures between sampling stations/ depth layers. Community structures of each sampling were visualized using non-metric dimensional scaling (NMDS) based on the Bray-Curtis dissimilarity index of the copy numbers of each OTUs. Similarity tests between two variables (sampling stations and depth layers) were conducted using nonparametric multivariate analysis of variance (NPMANOVA). All statistical analyses were performed using RStudio version 1.3.1093 25  For this study, we obtained 4,959,721 MiSeq reads, of which 3,597,424 were high-quality, merged reads. Number of non-standard fish Miseq reads were 2,170,008 out of the 3,597,424 (60.3%). In the MiSeq run, 18.5-98.6% of sequence reads out of the total reads were from non-standard fish DNAs of field samples, and 0% (= 0/159,453) was from those of a field negative control. The final list included 92 OTUs distributed across 49 families and 83 genera (Table S2).
Comparison of eDNA metabarcoding with fish catch record. The fish catch record of a set net indicated that 289 thousand kg of total fish biomass belonging to 41 classifications (38 species, one genus, and two families) was caught in May 2018. The catch biomass was highly variable between the classifications. For example, 255 thousand kg of Scomber japonicus was recorded, whereas Fistulariidae was 0.2 kg only (Table S3). Of the 41 classifications caught by the set net, eDNA metabarcoding detected 23 species. Three classifications of the 41 classifications were not recorded at the species level in the list of the fish catch of the set net, while the metabarcoding cannot assign OTUs at the species level for the Scomber japonicus and S. australasicus because they are so closely related. Thus, excluding these five species, our metabarcoding eventually detected 58.3% (= 21/36) of the fishes caught by the set net. Among the five dominant OTUs detected by eDNA metabarcoding described below (B. splendens, P. trilineatum, Scomber spp., P. major, and T. japonicus), four OTUs except B. splendens appeared in the list of the fish catch.
Relationship between the copy numbers and sequence reads of the standard DNA. The relationships between the sequence reads and copy numbers of standard DNAs were examined by linear regression analysis (Fig. 2). Within each sample, the sequence reads linearly and positively correlated with the copy numbers of standard DNAs (Fig. 2a). R 2 values of the regression lines ranged from 0.71 to 0.99, and more than 75% of the samples had values over 0.9 (Fig. 2b). These results suggested that the number of sequence reads was generally proportional to the number of copies of DNA initially added to the first PCR reaction within a single sample. The slopes of the linear regressions varied depending on the sample, suggesting that the amplification www.nature.com/scientificreports/ efficiency or removal amounts of DNAs during the purification of PCR products were dependent on the individual sample. These sample-specific slopes were used to convert sequence reads to the copy numbers of fish OTUs. Based on the estimated copy numbers of each OTU, five dominant OTUs were Beryx splendens, Parapristipoma trilineatum, Scomber spp. (S. japonicus or S. australasicus), Pargus major, and Trachurus japonicus (Fig. 1c-g). Levels of contaminations were assessed using the field negative control. DNA copy numbers detected in the field negative control (n = 1) were 0% (= 0/67,294) of mean DNA copy numbers of field-positive samples, suggesting that there was no amount of contamination during the library preparation process.
Site variation in fish eDNA quantities, the number of OTUs, and echo intensity. The estimated copy numbers of fish OTUs were summed to estimate the total, demersal, and pelagic fish eDNA quantities. The average eDNA copy numbers of total fish, demersal fish, and pelagic fish were 64.3, 44.6, and 19.7 copies/ ml water, respectively, while those of five dominant OTUs, B. splendens, P. trilineatum, Scomber spp., P. major, and T. japonicus, were 20.3, 13.8, 12.9, 5.5, and 3.2 copies/ml water, respectively. Likelihood ratio tests showed that the eDNA copy numbers significantly varied among sampling stations for total, demersal, pelagic fishes, P. trilineatum, Scomber spp., P. major, and T. japonicus, whereas no significant variation was detected for B. splendens (Table S4). Tukey-HSD tests indicated that significantly higher total and demersal fish eDNA copy numbers were observed at AR1 than at W150 and E 500-750, while the pelagic fish eDNA copy number was higher at AR1 than at E750 (Fig. S1). Also, the total, demersal, and pelagic fish eDNA quantities were significantly higher at AR2 than at E750. Among dominant fish OTUs, P. trilineatum and T. japonicus eDNA copies were significantly higher at AR1 than at almost all the surrounding stations, and P. major eDNA copies were higher at AR1 than at W150 and E 500-750. Scomber spp. eDNA quantities were higher at AR2 than at W150 and E 500-750. Echo intensity (S V ) ranged from 5.46 × 10 -9 to 5.20 × 10 -6 , and a likelihood ratio test indicated significant variation in echo intensity among the stations. Like the fish eDNA quantities, the echo intensity was higher at AR1 than at the surrounding stations. Meanwhile, the number of OTUs detected at each station ranged from 9 to 27, with a mean of 16 OTUs, and there was no significant variation among the stations (Table S4).

Effects of distance from an artificial reef and sampled depth layer on fish eDNA quantities and the number of OTUs.
Distance from the AR1 and depth layer were included in the Gamma component of the minimum BIC model for total fish, demersal fish, pelagic fish, and Scomber spp. eDNA quantities (Table 1). Estimated parameters indicated that these eDNA quantities decreased with increasing distance from the AR1 and the eDNA quantities were higher in the bottom samples than in the middle ones (Fig. 3). For Parapristipoma trilineatum, Pargus major, and Trachurus japonicus eDNA quantities, distance from the AR1, depth layer, and their interaction were included in the Gamma component of the minimum BIC model. Estimated parameters indicated that these eDNA quantities also decreased with increasing distance from AR1, but the intercepts and slopes of the relationship with distance differed between the layers. The eDNA quantities of P. major were higher in the bottom layer at 0 m from AR1, while those of P. trilineatum and T. japonicus were higher in the middle samples at 0 m from AR1 (Fig. 3). Slopes for P. trilineatum, and T. japonicus eDNA quantities were steeper in the middle layer, while that for P. major eDNA was steeper in the bottom layer and P. major eDNA decreased with distance from AR only in this layer. Only the depth layer was included in the Gamma component of the minimum BIC model for the number of OTUs and Beryx splendens eDNA quantity, and these values were higher in the bottom layer than in the middle layer (Table 1 and Fig. 3).
Fish eDNA quantities in relation to echo intensity and sampled depth layer. Echo intensity and sampled depth layer were selected as explanatory variables in the Gamma component of the minimum BIC model for total fish, demersal fish, pelagic fish, P. trilineatum, Scomber spp., and T. japonicus eDNA (Table S5). Estimated parameters indicated that these eDNA quantities increased with increasing echo intensity and the eDNA quantities were lower in the middle samples than in the bottom ones (Fig. 4). For P. trilineatum and P. major eDNA quantities, echo intensity, depth layer, and their interaction were included in the Gamma component of the minimum BIC model. Estimated parameters indicated that these eDNA quantities increased with increasing echo intensity, but the intercepts and slopes of the relationship with distance differed between the layers. The eDNA quantities of P. trilineatum were higher in the middle layer, while those of P. major were higher in the bottom samples. A slope for P. trilineatum was steeper in the middle layer, while that for P. major eDNA was steeper in the bottom layer. However, only depth layer was selected as an explanatory variable for B. splendens eDNA and its eDNA quantities were higher in the bottom layer. These minimum BIC models had high predictive power (R 2 ≥ 0.55) for these eDNA quantities ( Fig. 4 and Table S5).
Fish community structure. Community compositions of dominant fish OTUs at each site as well as each depth layer are shown in Fig. 5. P. trilineatum, P. major. B. splendens and T. Japonicus were dominant OTUs at AR1, while B. splendens and Scomber spp. were dominant at AR2. B. splendens and Scomber spp. were also dominant at surrounding stations such as E150, E500, W500, and W750. A comparison between the depth layers indicated that B. splendens, Scomber spp., and P. major were dominant in the bottom layer, while P. trilineatum was dominant in the middle layer. There was dissimilarity in the NMDS in fish community structure between the ARs and other surrounding stations (Fig. 6, NMDS stress = 0.096). PERMANOVA detected significant differences in community structure between the ARs (AR1 and AR2) and surrounding stations (E and W) as well as between the middle and bottom sampling layers (ARs vs surrounding stations: F = 3.041, p = 0.005; Middle vs bottom: F = 4.544, p = 0.001).

Discussion
Field studies have begun to use an eDNA approach for monitoring fish density or diversity in coastal ecosystems 6,7,9,10 , but, as far as we know, there have been no studies using this method to assess artificial reef (AR) effects on marine fishes. We found that eDNA metabarcoding efficiently detected the fish assemblages around ARs of Tateyama Bay. The eDNA metabarcoding detected 92 OTUs from 32 eDNA samples collected within 8 h on May 23, 2018 (Table S2), while 41 classifications (38 species, one genus, and two families) were recorded in the fish catch of the net set near the study site in May 2018 (Table S3). Of the 36 species recorded in the fish catch over a month, eDNA metabarcoding detected 58.3% of species (i.e., 21 species/36 species). This detection rate of metabarcoding relative to the conventional approach was within the range of previous eDNA metabarcoding studies with MiFish primers (45-70%) 6,33,34 . Also, 14 species of the recorded species in the fish catch were represented by < 10 kg over the one month. Meanwhile, our metabarcoding detected 77.3% (= 17 species/22 species) of fish species for which ≥ 10 kg were recorded in the fish catch (Table S3). Four dominant OTUs in terms of eDNA quantities (P. trilineatum, Scomber spp., P. major, and T. japonicus) were also a higher proportion of fish catch of the set net (Table S3). Meanwhile, the eDNA metabarcoding detected 68 fish OTUs that were not collected by the set net. This difference in the detected fish community between the eDNA metabarcoding and set net may be partly because the set net could not collect small-sized fish, and non-targeted fishes (e.g., Myctophidae and Gobiidae). These results indicated that our eDNA metabarcoding approach had a high and reasonable resolution for examining a fish community in the study site.
Results from eDNA metabarcoding analysis provide taxonomic table with the number of sequence reads, which are greatly affected by the amplification efficiency of the PCR primers and by PCR cycles 12 . Therefore, comparisons of the number of sequence reads among different samples and studies are not straightforward. However, we overcame this limitation by quantifing eDNA copy numbers of multiple fish species from their sequence reads adding internal standard DNAs at several levels of concentration (i.e., known amounts of short DNA fragments from fish species that have never been observed in a sampling area) to the eDNA samples 7 . We observed linear relationships between the number of sequence reads and the quantities of internal standard DNAs (Fig. 2). This result suggested that the sequence reads were proportional to the copy numbers of DNA within a single www.nature.com/scientificreports/ sample, indicating that the conversion of sequence reads using a sample-specific regression slope would provide an estimation of DNA copy numbers in the DNA extract (i.e., estimated DNA copy number = sequence reads/ sample-specific regression slope). This finding is consistent with previous studies showing that this procedure provides a reasonable estimate of copy numbers of multispecies eDNAs 7,16 . Therefore, this procedure enabled us to compare eDNA quantities between different samples and studies for multiple fish species. Meanwhile, the qMiSeq also has limitations, i.e., this method does not take account of the variation in amplification efficiency among species due to their differences in number of primer mismatches and barcode length. Although we found no differences in the primer mismatches among the nine dominant species of this study (Table S6) and consider there was no amplification bias at least in terms of primer compatibility among these species, we generally have to caution against comparisons of DNA copy numbers between different species using the qMiSeq. This study efficiently estimated the eDNA quantities of multiple fish species around an AR. The qMiSeq found that the eDNA quantities of total, demersal, pelagic fishes, and four dominant OTUs, except B. splendens, were higher at AR1 than adjacent sites. In addition, all eDNA quantities, except B. Splendens, sharply decreased with increasing distance from the AR. In particular, those of P. trilineatum, P. major, and T. japonicus clearly decreased even 150 m from AR1 (e.g., W150 and E150). These results are consistent with previous studies showing a higher density of fishes in ARs than reference sites 2,36,37 . In fact, four dominant OTUs, P. trilineatum, Scomber spp., P. major, and T. japonicus, were known to appear and aggregate in ARs 2,20,36,37 . Only B. splendens was not reported to be dependent on ARs, although it is highly targeted by fisheries and usually caught by vertical longline at midnight in a deeper zone of this region (depth > 200 m, personal communication, M. Sato). A review paper of the biology of B. splendens indicated that the fisheries grounds of this species are 100-500 m in depth around the mouth of Tokyo Bay, but their juveniles are sometimes collected in a shallower area (< 100 m) 38 . Although we cannot exclude the possibility that eDNA from adult B. splendens could be detected in this study, its higher eDNA quantities in the stations of 70-80 m depth (e.g., AR2, AR1 and E150) would be likely to reflect the presence of its juveniles. In future studies, combining eDNA analysis with other methods such as observation by Depth layer Figure 3. The relationships between distance from AR1 as well as number of OTUs and eDNA copies of total fish, demersal fish, pelagic fish, Beryx splendens, Parapristipoma trilineatum, Scomber spp., Pagrus major, and Trachurus japonicus. The red and turquoise blue lines indicate median values in the bottom and middle layers respectively, while each color circle indicates raw values in the respective layers. The gray shaded area represents a 95% confidence interval, as determined from the generalized linear model. We expected the number of OTUs would be higher at ARs than the surrounding stations but did not observe such a pattern. In the following survey in June 2019, we found a larger number of fish species around AR1 than surrounding stations in a GoPro video (personal observation, N. Inoue, M. Sato, N, Furuichi), which contradicted our eDNA results. Because MiSeq sequencing produces a limited amount of reads and this limited number of reads were distributed among the OTUs in the flowcell, we suspect that higher eDNA concentrations of dominant OTUs at ARs consumed a large portion of the sequence reads, and species with low eDNA concentrations were not detected by the metabarcoding in this study. On the contrary, the qMiSeq successfully converted sequence reads into DNA copy numbers for multiple species using a sample-specific regression slope (Fig. 2). Ushio et al. 7 already showed that the numbers of eDNA copies estimated by qMiSeq and qPCR were significantly and positively correlated with each other for total fish and dominant species. Therefore, we consider that a limited number of reads in MiSeq sequencing do not prevent from spatial comparisons of eDNA copies among sampling locations especially for total fish and dominant species.
The eDNA metabarcoding also showed different fish species composition between the AR and surrounding stations based on NMDS and PERMANOVA. The eDNA quantities of fishes, including dominant OTUs (pelagic fishes: Scomber spp., Etrumeus teres, Sardinops melanostictus; demersal fishes: B. splendens, P. trilineatum, P. major, Sardinops melanostictus, Sacura margaritacea, Thamnaconus modestus) were more abundant in ARs, while those of some demersal fishes such as flathead (Platycephalus sp.2), Asian sheephead wrasse (Semicossyphus reticulatus), barred knifejaw (Oplegnathus fasciatus), and seabream (Dentex tumifrons) were higher in the surrounding  www.nature.com/scientificreports/ stations. Different fish responses to AR structures may result in differences in the fish community compositions between the ARs and surrounding stations. Previous studies have found positive relationships between the field density of target species and eDNA concentrations in semi-closed marine systems like estuaries 9,14 ; our results found such associations in an open marine system. We found eDNA quantities of total, demersal, pelagic fishes, and four dominant species, except B. splendens, positively correlated with the echo intensity measured by echo sounder. Our study site faces toward the open sea (i.e., Pacific Ocean), and the flow was stronger than in a previous study (5.0-25.0 cm/s in our study and 0.4-13.5 cm/s in a semi-closed bay at Maizuru, Japan) 39 . Therefore, the flow is likely to transport more eDNAs from sampling stations in this study, but we could also obtain a strongly localized eDNA signal at sampling stations. One exception is the eDNA quantity of B. splendens, which did not positively correlate with echo intensity in both the middle and bottom layers. This finding was probably because small-sized juveniles of this species may be detected by eDNA survey but not by acoustic survey. Overall, our result indicates that eDNA quantities reflect spatial distributions of marine fish in an open coastal environment.
We found the eDNA quantities and species number were higher in the bottom sampling layers than in the middle layers. Multivariate analyses, including NMDS and PERMANOVA, also showed different fish compositions between the bottom and middle layers. Previous eDNA studies have also indicated variations in eDNA quantities, the number of species, and community compositions between different sampling depths 6,9,40,41 . Some of these studies using eDNA methods suggested that water column stratification and the different species' responses to it cause the observed different composition between sampling depths 40,41 . Murakami et al. 39 also discussed that the vertical distribution of eDNA depends on both the depth of the source organisms and the behavior of eDNA, with the latter dependent on its chemical condition and physical factors in the environment. Although water column stratification between the middle and bottom layers was not observed in our study (Fig. S2), we consider that the preference of detected fish species for bottom habitats or the behavior of eDNA contributed to the higher eDNA quantities and OTU number in the bottom layer in this study.
Overall, our results demonstrate the effectiveness of qMiSeq to assess AR effects on multiple fish species. We found higher quantities of fish eDNA at ARs than at surrounding stations and positive correlations between these quantities and echo intensity, which indicates a highly localized eDNA signal in an oceanographically dynamic environment. This spatial resolution indicates that qMiSeq can improve techniques to monitor the distribution of multiple fish species. Another method, called 'quantitative sequencing (qSeq)' has also been developed for simultaneous quantifications of eDNA concentrations of multiple species 42 . This method adopts a stochastic labeling approach with random-tag sequences and can also quantify eDNA concentrations of fish species using the MiFish primer set 43 . The quantitative metabarcoding methods (i.e., qMiSeq and qSeq) can complement conventional methods, especially in an environment with areas where the water is too deep to adopt underwater visual census (> 40 m depth) or in no-take marine reserves where exploiting fish is prohibited.

Data availability
DDBJ accession numbers of the DNA sequences analyzed in the present study are DRA011864 (Submission ID), PRJDB11389 (BioProject ID) and SAMD00298249-SAMD00298281 (BioSample ID).