Quantitative assessment of multiple fish species around artificial reefs using environmental DNA metabarcoding

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) 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 quantitative MiSeq sequencing and echo sounder survey. Using the quantitative metabarcoding method, we found higher quantities of fish eDNAs at the ARs than at surrounding stations and different fish species compositions between them. Comparisons with eco sounder survey also showed positive correlations between fish eDNA concentration and echo intensity, which indicates a highly localized signal of eDNA at each sampling station. These results suggest that quantitative MiSeq sequencing is a promising technique to complement conventional methods to monitor distributions of multiple fish species.


8
Artificial reefs (ARs) have been deployed worldwide to create fisheries grounds. Since the early 3 9 1970s, restrictions on overseas fishing have influenced high seas fisheries outside the exclusive 4 0 economic zones (EEZ), and this situation led the Japanese government to implement a new 4 1 fishery policy to develop and maintain coastal fisheries. In 1971, the Japanese government started 1 8 2 initial base call calibrations on the MiSeq platform. The thermal cycle profile after an initial 3 min 1 8 3 denaturation at 95˚C was as follows (35 cycles): denaturation at 98˚C for 20 s; annealing at 65˚C 1 8 4 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 1 8 6 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 SPRIselect (Beckman Coulter inc.) to remove dimers and monomers following the 1 9 0 manufacturer's protocol. The second-round PCR (2nd PCR) was carried out with a 24 µl reaction volume containing 12 1 9 2 µl of 2 × KAPA HiFi HotStart ReadyMix, 2.8 µl of each primer (5 µM), 4.4 µl of sterilized 1 9 3 distilled H2O, and 2.0 µl of template. We used the following two primers to append the 1 9 4 dual-index sequences (8 nucleotides indicated by Xs) and flowcell-binding sites for the MiSeq profile after an initial 3 min denaturation at 95˚C was as follows (12 cycles): denaturation at 98˚C 2 0 0 for 20 s; combined annealing and extension at 72˚C for 15 s, with a final extension at 72˚C for 5 2 0 1 min. The concentration of each second PCR product was measured by quantitative PCR using TB 2 0 2 Green Fast qPCR Mix (Takara inc.). Each sample was diluted to a fixed concentration and 2 0 3 combined (i.e., one pooled 2nd PCR product that included all samples). The pooled 2nd PCR 2 0 4 product was size-selected to approximately 370 bp using BluePippin (SAGE SCIENCE). The 2 0 5 size-selected library was purified using the Agencourt AMPure XP beads, adjusted to 4 nM by 2 0 6 quantitative PCR using TB Green Fast qPCR Mix (Takara inc.), and sequenced on the MiSeq 2 0 7 platform using a MiSeq v2 Reagent Kit (2 × 150 bp). 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 2 1 2 implemented in Claident 17 . We adopted this process rather than using FASTQ files demultiplexed 2 1 3 by the Illumina MiSeq default program in order to remove sequences whose 8 mer index 2 1 4 positions included nucleotides with low-quality scores (i.e., Q-score < 30).

1 5
The processed reads were subjected to a BLASTN search against the full NCBI database. If the 2 1 6 sequence similarity between queries and the top BLASTN hit was < 98.5% and the sequence 2 1 7 7 length was less than ≤ 150 bp, the unique sequence was not subjected to the following analyses.

1 8
After BLASTN searches, assembled sequences assigned to the same species were clustered, and 2 1 9 we considered the clustered sequences as operational taxonomic units (OTUs). To remove 2 2 0 possible contaminants, we removed the OTUs whose sequence reads were < 0.05 % of the total 2 2 1 reads as a noise for each sample. To confirm the feasibility of MiFish metabarcoding, we compared the fish community of MiFish 2 2 5 metabarcoding with a catch record of a net set near the study site (Fig.1). The detection rate was 2 2 6 defined as the proportion of species detected by metabarcoding compared to the species collected 2 2 7 by the net. Because our eDNA survey was conducted in May 2018, we used catch-record data 2 2 8 from the same period (May 2018). According to Ushio et al. 7 and Ushio 15 , the procedure to estimate DNA copy numbers consisted of sequence reads of non-standard fish DNAs to estimate the copy numbers using the result of the  Linear regressions were used to examine how many sequence reads were generated from one 2 3 7 fish DNA copy through the library preparation process. Note that a linear regression analysis 2 3 8 between sequence reads and standard DNAs was performed for each sample, and the intercept and thus thirty-three regression slopes were estimated in total.

4 3
The sequence reads of non-standard fish DNAs were converted to copy numbers using 2 4 4 sample-specific regression slopes estimated by the above regression analysis. The number of 2 4 5 non-standard fish DNA copies was estimated by dividing the number of MiSeq sequence reads by 2 4 6 a sample-specific regression slope (i.e., the number of DNA copies = MiSeq sequence 2 4 7 reads/regression slope). A previous study demonstrated that these procedures provide a 2 4 8 reasonable estimation of DNA copy numbers using high-throughput sequencing 7 . We estimated relative fish density using a quantitative echo sounder following the acoustic survey  For this study, we obtained 4,862,616 MiSeq reads, of which 3,528,543 were high-quality, 3 1 2 merged reads. Number of non-standard fish Miseq reads were 2,153,777 out of the 3,528,543 3 1 3 (61.0%). In the MiSeq run, 18.8%-98.6% of sequence reads were from non-standard fish DNAs 3 1 4 of field samples, and 0% was from those of a field negative control. The final list included 95 3 1 5 OTUs distributed across 49 families and 83 genera (Table S1). The fish catch record of a set net indicated that 289 thousand kg of total fish biomass belonging to OTUs except B. splendens appeared in the list of the fish catch. The relationships between the sequence reads and copy numbers of standard DNAs were  Levels of contaminations were assessed using the field negative control. DNA copy numbers 3 4 5 detected in the field negative control (n = 1) were 0% of mean DNA copy numbers of  Site and vertical 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 3 5 1 pelagic fish eDNA quantities. The average eDNA copy numbers of total fish, demersal fish, and 3 5 2 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, 3 5 4 12.9, 5.5, and 3.2 copies/ml water, respectively. Likelihood ratio tests showed that the eDNA 3 5 5 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 3 5 7 detected for B. splendens (Table S3). Tukey-HSD tests indicated that significantly higher total 3 5 8 and demersal fish eDNA copy numbers were observed at AR1 than at W150-750 and E 500-750, while the pelagic fish eDNA copy number was higher at AR1 than at E750 (Fig. S1). Also, the 3 6 0 total, demersal, and pelagic fish eDNA quantities were significantly higher at AR2 than at E750. For dominant fish OTUs, P. trilineatum, P. major, and T. japonicus eDNA copies were 3 6 2 significantly higher at AR1 than at other stations, while Scomber spp. eDNA quantities were 3 6 3 higher at AR1 and AR2 than at E750. Echo intensity (S V ) ranged from 4.53×10 -8 to 1.38×10 -6 , and Field studies have begun to use an eDNA approach for monitoring fish density or diversity in 4 0 4 coastal ecosystems 6,7,9,10 , but, as far as we know, there have been no studies using this method to 4 0 5 assess artificial reef (AR) effects on marine fishes. We found that eDNA metabarcoding japonicus) were also a higher proportion of fish catch of the set net (Table S1). Meanwhile, partly because the set net could not collect small-sized fish, and non-targeted fishes (e.g.,  We observed linear relationships between the number of sequence reads and the quantities of Miseq also has limitations; i.e., this method does not take account of the variation in amplification 4 2 9 efficiency among species. Therefore, we have to caution against comparisons of DNA copy  This study efficiently estimated the eDNA quantities of multiple fish species around an AR.

3 2
Quantitative MiSeq sequencing found that the eDNA quantities of total, demersal, pelagic fishes, showing a higher density of fishes in ARs than reference sites 2,30,31 . In fact, four dominant OTUs, 4 3 8 P. trilineatum, Scomber spp., P. major, and T. japonicus, were known to appear and aggregate in 4 3 9 ARs 2,31-33 . Only B. splendens was not reported to be dependent on ARs, although it is highly We expected the number of OTUs would be higher at ARs than the surrounding stations, but 4 5 0 did not observe such a pattern. In the following year's survey, we found a larger number of fishes  MiFish metabarcoding also showed different fish species composition between the AR and    simultaneous detection of multiple fish species from environmental DNA and other samples.  benchmark using barcode sequences of bacteria, Archaea, animals, fungi, and land plants.