The field survey was performed in Tateyama Bay, central Japan, in the proximity of the Kuroshio warm current facing the Pacific Ocean (Fig. 1). Many ARs have been implemented in this area. In our study, we focused on the highest steel AR (UT-304, Shinko Kenzai Ltd., Tokyo, Japan) at 34°59'44"N, 139°47'30"E, which is 30 m in height deployed at a depth of 75 m in the waters of Tateyama Bay and has a higher density of fish species than other ARs (Sato et al. 2021). Water sampling for eDNA analysis was undertaken on June 3 (Day 1) and 4 (Day 2), 2019, on board the R/V Taka-maru (Japan Fisheries Research and Education Agency: FRA). Sampling stations were set up at the AR and at three distant points (0.1, 0.2, and 0.4 nautical miles; 0.1 miles ≈ 185 m) extending from the AR in four directions (Fig. 1). Sampling was conducted on Days 1 and 2 at a total of 13 stations (due to the lack of survey time by the vessel, excluding the southern three stations on Day 2). Details of the sampling procedure are shown in Suppl. material 1: Table S1.

Figure 1.

Maps showing the wider sampling area encompassing Japan (A), Tokyo Bay (B), and sampling stations (C). Schematic diagram of the water sampling method (D).

At each sampling station, 10 L of seawater was collected from a depth of 40 m, which is about the same depth as the AR except for St. S3, using one cast of two Niskin water samplers (5 L × 2 samples) (Fig. 1). At St. S3, whose depth is shallower than 40 m, seawater was collected from 20 m depth using the same method with Niskin water samplers. At the AR, we additionally collected 10 L of seawater from the surface using a bucket and 10 L from the bottom (70 m depth) using one cast of two Niskin samplers of the same type. Two 2-L samples were subsampled from the 10-L seawater sample using a measuring bottle, and the remaining 6 L of seawater was used to pre-wash the devices used. Then, in a laboratory on the R/V Taka-maru, following the method described by Miya et al. (2016), two 2-L samples were immediately filtered using a combination of Sterivex filter cartridges (nominal pore size = 0.45 μm; SVHV010RS, Merck Millipore, Darmstadt, Germany) through an aspirator. After filtration, 1.5 ml of RNAlater (ThermoFisher Scientific, Waltham, MA, USA) was injected into the cartridge to prevent eDNA degradation. The Niskin water samplers and bucket were bleached before each water collection using a commercial bleach solution, while filtering devices were bleached after each filtration. Two L of MilliQ water were filtered with a filter funnel and measuring cup for each day as a field negative control to test for possible contamination. The filter cartridges were placed in a freezer and stored at below –30 °C immediately after filtration until eDNA extraction.

The eDNA was extracted and purified following the method developed by Miya et al. (2016). After removing RNA later from inside the cartridge using a centrifuge, a mixture 60 μl of proteinase-K, 660 μl of phosphate buffered saline, and 600 μl of buffer AL was injected into the cartridge and incubated 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, Venlo, Netherlands) following the manufacturer’s protocol. After purification, DNA was eluted in 100 ul AE buffer and were frozen at –20 °C.

For quantitative MiSeq sequencing, five artificially designed internal standard DNAs (Ushio et al. 2018; Ushio 2019) were used to calculate standard curves to estimate the DNA copy numbers. The copy numbers of the five standard which were artificially designed and synthetic, A to D DNAs were adjusted to 100, 50, 25, 12.5, and 2.5 copies μl -1. Sequences of the five standard DNAs are shown in Suppl. material 1: Table S2. Two PCR‐level negative controls (i.e. with and without internal standard DNAs) were employed for the MiSeq run to monitor contamination during the experiments.

The first-round PCR (1^st PCR) was carried out with a 12-μl reaction volume containing 6.0 μl of 2× KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, WA, USA), 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. The standard DNA mix was included for each sample. The final concentration of each primer was 0.3 μM. A mixture of the following four PCR primers modified from original MiFish primers was used (Miya et al. 2015): MiFish-U-forward (5´-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT NNN NN G TCG GTA AAA CTC GTG CCA GC) 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´). The thermal cycle profile after an initial 3-min denaturation at 95 °C was as follows (35 cycles): denaturation at 98 °C for 20 sec; annealing at 65 °C for 15 sec; and extension at 72 °C for 15 sec, with a final extension at the same temperature for 5 min. Eight replications were performed for the 1^st PCR, and products were pooled in a single 1.5-ml tube. The pooled products were purified and size-selected for 200–400 bp using a SPRIselect (Beckman Coulter, Fullerton, CA, USA) to remove dimers and monomers following the manufacturer’s protocol.

The second-round PCR (2^nd 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 H2O, and 2.0 μl of template. We used the following two primers to append the dual-index sequences (eight nucleotides indicated by Xs) and flowcell-binding sites for the MiSeq platform (5´ ends of the sequences before the eight Xs): 2^nd-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 2^nd-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 2^nd-PCR product was measured by quantitative PCR using TB Green Fast qPCR Mix (TaKaRa, Otsu, Japan). Each sample was diluted to a fixed concentration and combined (i.e. one pooled 2^nd-PCR product that included all samples). The pooled 2^nd-PCR product was size-selected to approximately 370 bp using BluePippin (Sage Science, Beverley, MA, USA). The size-selected library was purified using Agencourt AMPure XP beads, adjusted to 4 nM by quantitative PCR using TB Green Fast qPCR Mix (TaKaRa, Otsu, Japan), and sequenced on the MiSeq platform using a MiSeq v2 Reagent Kit (Illumina, San Diego, CA, USA).

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 and the processed reads were subjected to a BLASTN search against the full NCBI database using Claident (Tanabe and Toju 2013). Reads where the sequence similarity between queries and the top BLASTN hit was < 98.5% and the sequence length was less than ≤ 150 bp were not used for the subsequent analyses. After BLASTN searches, the assembled sequences assigned to the same species were clustered, and the clustered sequences were treated as operational taxonomic units (OTUs). The OTUs whose sequence reads were < 0.05% of the total reads were considered a noise for each sample and thus excluded to remove possible contaminants. In the present study, after OTUs read by the above-mentioned methods were removed, the remaining OTUs were defined as a high-quality read. The numbers of DNA copies of each sample were estimated using methods described by Ushio et al. (2018) and Ushio (2019): first, the linear regression of five internal standard DNA copies and read numbers were used to estimate the standard curve, and then the curve was used to convert the read of each OTU to DNA copy numbers. At least three out of five internal standard DNA copies were detected samples that were used for subsequent analyses in the present study.

The sea current field, i.e. the east-west velocity (u, eastward positive) and the north-south velocity (v, northward positive), was measured at multiple layers using a ship-mounted 300-kHz ADCP (Teledyne RD Instruments, Poway, CA, USA). In subsequent analyses, velocity data from the layer at 47 m depth (excluding St. S3 where the depth was 23 m) were used. Both layers were close to the middle layer sampling depth zone. Velocity data were time-averaged for 7 min, i.e. from 3 min before the 1-min sampling time to 3 min after the sampling time (Suppl. material 1: Table S1). The parameter Dir (index of current direction relative to the AR) was introduced to evaluate the effect of the sea current direction relative to the AR on the distribution of eDNA density. The index Dir quantifies to what extent the current at each sampling station is directed toward the AR.

As illustrated in Fig. 2, the index Dir was defined as

$Dir\equiv \overrightarrow{r}·\overrightarrow{U}/U$ (1)

where $\overrightarrow{r}$ is the unit vector indicating the direction from the AR to each sampling station, $\overrightarrow{U}$ = (u, v) = (U cosθ, U sinθ) is the velocity vector, U = (u² + v²)^0.5 is the absolute value of velocity, θ is the current direction estimated from the time-averaged u and v data (0 degrees for the eastward direction and counter-clockwise positive), and the operator “.” denotes the inner product. By incorporating $\overrightarrow{r}$ = (0, 1) for the northern stations, (0, -1) for the southern stations, (1, 0) for the eastern stations, and (-1, 0) for the western stations into equation (1), respectively, index Dir is expressed more specifically as sin (θ) (northern stations), –sin (θ) (southern), cos (θ) (eastern), and -cos (θ) (western). The value of Dir ranges from -1.0 to 1.0, where negative values (blue in Fig. 2) indicate upstream currents relative to the AR and positive values (red in Fig. 2) indicate downstream currents.

Figure 2.

Schematic diagram for the calculation of the index Dir. Dir is defined as the inner product of the unit vector indicating the direction from the AR to each sampling statio ($\overrightarrow{r}$) and the velocity vecto ($\overrightarrow{U}$) . By assuming $\overrightarrow{r}$ = (0, 1) for the northern stations, (0, -1) for the southern stations, (1, 0) for the eastern stations, and (-1,0) for the western stations, Dir is expressed more specifically as sin (θ) (northern), -sin (θ) (southern), cos (θ) (eastern), and -cos (θ) (western). The value of Dir ranges from -1.0 to 1.0, where the negative (blue color) and positive (red color) values respectively indicate upstream and downstream directions from the AR.

As described in Results, the top four species with the largest amount of eDNA distribution were Japanese anchovy (Engraulis japonicus), threeline grunt (Parapristipoma trilineatum), chub mackerel (Scomber spp), and horse mackerel (T. japonicus) (Table 1). To examine the relationship between the eDNA density distribution of the four dominant species, and their distance from the AR and the current fields around the AR, general linear modeling (GLM) was used to create current-distribution models. These models, which assumed a normal distribution with an identity link function with a log-based eDNA density (log (eDNA+0.1)) as the response variable, and day (Day), distance from the AR (Dist), index of the relative current direction to the AR (Dir), velocity (Vel), and interaction between Dir and Vel (Dir: Vel) as explanatory variables, were created for the initial models. Akaike’s information criterion (AIC) and Bayesian information criterion (BIC) were then used to select the best model among the initial models. The Wald test was used to determine the significance for coefficients of the explanatory variables of the best models. The likelihood ratio test was used to compare the the best and null models. All statistical analyses were performed using R (version 4.0.4; R Development Core Team 2020).

The MiSeq paired‐end sequencing of the 58 libraries, which comprised 54 field samples (Day 1: 30 samples, Day 2: 24 samples) and both a negative field sample and a PCR-level control sample on each day for this study, yielded a total of 7,054,162 reads with 94.0% base calls containing Phred quality scores of over 30.0 (Q30; error rate = 0.1%, or base call accuracy = 99.9%), and 6,877,676 (97.5%) were high-quality reads of which 3,735,725 (54.3%) were non-standard fish sequences, and 135 operational taxonomic units (OTUs), which included at least 117 genera, were detected (Suppl. material 1: Table S3). Excluding the standard DNA read numbers, the sequence reads of the two negative field controls were 0 and 1143 (0% and 7.1% of the mean DNA copy numbers of field-positive samples), and the sequence reads of the two PCR-level negative controls were 0 and 391 (0% and 0.57% of the mean DNA copy numbers of field-positive samples). This result indicated that serious contamination during field sampling and/or laboratory experiments did not occur or was insignificant. Data sets of three samples – one of which was collected from St. W2, one from the water surface at the AR, and one from its middle layer – were excluded from subsequent analyses, because three out of five standard DNAs were not detected in these samples (i.e. the sequence reads of these samples could not be converted to accurate eDNA densities). The sample-specific relationship between the sequence reads and the copy numbers of standard DNAs were examined using a linear regression. The mean r² value of the regression lines was 0.96, and more than 84% of the samples had values greater than 0.9 (Fig. 3). This indicates that the read numbers of each OTU in each sample could be converted to a DNA copy number with minimal error using the method described by Ushio et al (2018).

Figure 3.

Relationships between sequence reads and copy numbers of standard DNAs according to Ushio et al. (2018). Relationships between sequence reads and copy numbers of standard DNAs with maximum (open squares), median (black triangles), and minimum (open circles) slopes (A). Distribution of r² values of the regression lines (B).

The top ten species with the eDNA distribution are shown in Table 1. The present study focused on four dominant species that accounted for a total of 84.7% of the total number of DNA copies: 74.2% was for Japanese anchovy (E. japonicus); 5.7% for threeline grunt (P. trilineatum); 2.6% for chub mackerel (Scomber spp), which could not be identified as either S. japonicus or S. australasicus from the nucleotide sequence; and 2.2% for horse mackerel (T. japonicus) (Table 1).

Based on the ADCP data, current velocity and the direction at each sampling station are shown in Fig. 4. Current velocity on Day 1 was larger than that on Day 2, and northwest and northeast currents were prevalent, while on Day 2, west and southwest currents were prevalent. On Day 1, the indices of current direction relative to the AR (Dir) were negative (Dir < 0) at the three southern stations (S1, S2 and S3), one eastern station (E3), and two western stations (W2 and W3). Waters to the south and west of the AR were thus upstream of the current on Day 1. On Day 2, Dir was negative (Dir < 0) at one northern station (N2) and the three eastern stations (E1, E2 and E3). Waters to the east of the AR were thus upstream of the current on Day 2.

Figure 4.

Current direction and velocity at each station on Day 1 (A) and Day 2 (B). Blue and red arrows indicate the mean of 7-min ADCP data and show the upstream (Dir < 0) and downstream (Dir > 0) currents, respectively. Gray arrows show the current and velocity mean of each minute.

Spatial distributions of the eDNA densities for the four species are shown in Suppl. material 2: Figs S1–S4. In the case of E. japonicus (Suppl. material 2: Fig. S1), eDNA density at the AR was 5 to 30 times higher than at the other stations except for two western stations (W2 and W3), which were upstream of the AR on Day 1 (Fig. 4). On Day 2, eDNA density was low overall, just 1/30 of the mean density on Day 1 at the AR. However, the eastern stations, which were upstream of the AR, had the highest eDNA densities. In the case of P. trilineatum (Suppl. material 2: Fig. S2), eDNA density on Day 1 was the highest at western St. W2, which was upstream of the AR compared to the other stations. On Day 2, eDNA densities for P. trilineatum as well as the other species at two eastern stations (E1 and E3) were higher than at the other stations. In the case of Scomber spp (S. japonicus or S. australasicus) (Suppl. material 2: Fig. S3) as well as E. japonicus and P. trilineatum, on Day 1, eDNA densities at two western stations (W2 and W3) were higher than at the other stations. On Day 2, eDNA densities for all four species at two eastern stations (E1 and E3) were higher than at the other stations except for Scomber spp at St. W1. In the case of T. japonicus (Suppl. material 2: Fig. S4), eDNA density on Day 1 at St. S1, which was upstream of the AR, was higher than it was at the other stations, and eDNA densities at the southern stations decreased with increasing distance from the AR.

To examine the relationship between eDNA density and current fields for the four species, we created current-distribution models based on the GLMs. The results of model selection and statistical parameters of best models are shown in Tables 2, 3, respectively.

In the case of E. japonicus, out of five explanatory variables in the first model (AIC=147.567; BIC=159.565), Day, Dist and Dir were selected as explanatory variables in the best models (AIC=145.076; BIC=153.644) by AIC selection. Furthermore, Dist was removed from the best models by BIC selection, and Day and Dir were selected as explanatory variables in the BIC-based best models (AIC=146.184; BIC=153.038). Best models selected by AIC and BIC were significantly different from the null models (likelihood ratio test, AIC model: DF (37, 40), -2logL=-32.632, p<0.001; BIC model: DF (38, 40), -2logL= -27.535, p<0.001). For the AIC-based best model, coefficients of two explanatory variables (Day and Dir) were significant (Wald test, Day: t=-3.120, p=0.004; Dir: t=-2.697, p=0.012, respectively) and Dist was not significant (t=-1.707, p=0.096). However, the deviance explained (%) of the AIC-selected model was 33.5%, which was higher than that of the BIC-selected model (28.3%), while the difference between the BIC in both models (delta BIC) was as small as 0.61, indicating an insignificant difference between both best models (likelihood ratio test, df (38,37), -2logL= 5.097, p>0.05). The AIC-based best model was thus used to predict eDNA density, as described below.

For both P. trilineatum and Scomber spp, out of five explanatory variables in the first models (P. trilineatum: AIC=179.985, BIC=191.980; Scomber spp: AIC=205.562, BIC=217.557), Dir was selected as an explanatory variable in the best models (P. trilineatum: AIC=173.644, BIC=178.785; Scomber spp: AIC=199.822, BIC=204.962) by both model selections.

In the case of T. japonicus, out of five explanatory variables in the first model (AIC=185.173; BIC=197.168), Dist and Dir were selected as explanatory variables in the best models (AIC=182.465; BIC=189.319) by AIC selection. Furthermore, Dist was removed from the best models by BIC selection, and Dir was selected as an explanatory variable in the BIC-based best models (AIC=183.528; BIC=188.669). The best models selected by AIC and BIC were significantly different from the null models (likelihood ratio test, AIC model: DF (38, 40), -2logL=-37.409, p<0.05; BIC model: DF (39, 40), -2logL= -24.286, p<0.05). For the AIC-based best model, the coefficient of Dir was significant (Wald test, t=-2.626, p=0.012) and that of Dist was not significant (Wald test, t=-1.717, p=0.094). However, the deviance explained (%) of the AIC-selected model in the case of T. japonicus was 18.1%, which was also higher than that of the BIC-selected model (11.1%), and the difference between BIC in both models (delta BIC), at 0.65, was almost equally small, again indicating an insignificant difference between both best models (likelihood ratio test, df (39,38), -2logL=13.123, p>0.05). As with E. japonicus, the AIC-based best model was used to predict eDNA density for T. japonicus, too.

Using the best models, eDNA density in relation to Dir was predicted for the four species as follows (Figs 5–8): in the case of E. japonicus (Fig. 5), we predicted the eDNA density for LSmean of Day when Dir varied from -1 to 1 at Dist 0.1, 0.2, and 0.4 miles; for P. trilineatum (Fig. 6) and Scomber spp (Fig. 7), we predicted the eDNA density when Dir varied from -1 to 1; and in the case of T. japonicus, we predicted the eDNA density when Dir varied from -1 to 1 at Dist 0.1, 0.2, and 0.4 miles (Fig. 8). A similar trend was observed for all species in the relationship between eDNA density and Dir, i.e. the eDNA density of all species reached the maximum Dir at -1.0, which refers to the upstream side, and eDNA density rapidly decreased with an increase of Dir. In addition to Dir, Day and Dist were selected for the AIC-based best model in the case of E. japonicus, and Dist was selected for the AIC-based best model in the case of T. japonicus (Table 2). The pattern of eDNA density distribution was similar for both species: the more upstream and the closer to the AR the distribution was, the denser the eDNA was predicted to be (Figs 5, 8).

Figure 5.

Relationship between the index of relative current direction to the AR (Dir) and eDNA copies (A) and predicted eDNA copies (B) of Engraulis japonicus. Dir values range from -1.0 to 1.0. Negative values (Dir < 0) indicate upstream currents and positive values (Dir > 0) indicate downstream currents. Open triangles, black circles, and gray squares of (A) show the number of eDNA copies at Dist 0.1, 0.2, and 0.4, respectively. The dotted, solid, and gray lines of (B) show the predicted number of eDNA copies at Dist 0.1, 0.2, and 0.4, respectively.

Figure 6.

Relationship between the index of relative current direction to the AR (Dir) and number of eDNA copies (A) and predicted number of eDNA copies (B) of Parapristipoma trilineatum. The solid line and dotted lines of (B) show the predicted number of eDNA copies and 95% confidence intervals, respectively.

Figure 7.

Relationship between the index of relative current direction to the AR (Dir) and number of eDNA copies (A) and number of predicted eDNA copies (B) of Scomber spp (S. japonicus or S. australasicus). The solid line and dotted lines of (B) show the predicted number of eDNA copies and 95% confidence intervals, respectively.

Figure 8.

Relationship between the index of relative current direction to the AR (Dir) and number of eDNA copies (A) and predicted number of eDNA copies (B) of Trachurus japonicus. Open triangles, black circles, and gray squares of (A) show the number of eDNA copies at Dist 0.1, 0.2, and 0.4, respectively. The dotted, solid, and gray lines of (B) show the predicted number of eDNA copies at Dist 0.1, 0.2, and 0.4, respectively.