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

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¹⁴. 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 developed by Miya et al.¹⁵, 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 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 standard DNAs) were employed for MiSeq run to monitor contamination during the experiments. The first-round PCR (1st PCR) was carried out with a 12-µl reaction volume containing 6.0 µl of 2 × KAPA HiFi HotStart ReadyMix (Roche, Basel, Switzerland), 0.7 µl of each primer (5 µM), 2.6 µl of sterilized distilled H₂O, 1.0 µl of standard DNA mix and 1.0 µl of template. Note that the standard DNA mix was included for each sample. The final concentration of each primer was 0.3 µM. We used a mixture of the following four PCR primers modified from original MiFish primers¹⁶: MiFish-U-forward (5′-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT NNN NNG TCG GTA AAA CTC GTG CCA GC-3′) and MiFish-U-reverse (5′-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC TNN NNN CAT AGT GGG GTA TCT AAT CCC AGT TTG-3′), MiFish-E-forward (5′-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT NNN NNG TTG GTA AAT CTC GTG CCA GC-3′), and MiFish-E-reverse (5′-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC TNN NNN CAT AGT GGG GTA TCT AAT CCT AGT TTG-3′). 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 (^circ)C was as follows (35 cycles): denaturation at 98 (^circ)C for 20 s; annealing at 65 (^circ)C for 15 s; and extension at 72 (^circ)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₂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 (^circ)C was as follows (12 cycles): denaturation at 98 (^circ)C for 20 s; combined annealing and extension at 72 (^circ)C for 15 s, with a final extension at 72 (^circ)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¹⁷. 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).

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⁶. Following Yamamoto et al.⁶, 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.⁷ and Ushio¹⁵, 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⁷.

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.⁹. 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¹⁸ 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 AR1, the second and third analyses were separately performed. Based on Fishbase²¹, each fish OTU was classified into demersal and pelagic species. For the analyses, we used the Poisson distribution with a log-link function for OTU number because it was count data, while a gamma distribution with a log-link was used for eDNA copy numbers of total fish, pelagic fish, demersal fish, and Scomber spp. (one dominant OTU) because these values are continuous and greater than zero. A mixture of binomial and gamma distributions, which is often called a “gamma hurdle model,” was also used as in previous studies^22,23 for eDNA copy numbers of remaining dominant OTUs because these values are continuous and contain zero values. Bayesian information criterion (BIC) was used to select the correct model among candidate models for the second and third analyses²⁴.

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²⁵ (R version 3.6.1²⁶) with several packages. The packages used were: (a) glmmTMB²⁷ for regression models; (b) multcomp²⁸ for Tukey-HSD test; (c) MuMIn²⁹ for BIC model selection, (d) vegan³⁰ for multivariate analysis; (e) ggplot2³¹ and (f) ggpubr³² for figure creation.

Source: Ecology - nature.com