The trawl survey was conducted from June 15 to 20, 2022. The survey sea area was around the Yellow River estuary, specifically including 15 stations between E119.100 and E119.686 and N37.500 and N38.220 ( Figure 1 and Table 1 ). A single trawler with a power of 45.78 kW was used in the survey. The bottom trawl net port had a height of 3.0 m, a width of 10 m, and a circumference of 800 m, and the mesh size of the cod end was 20 mm. The average trawling speed was 2 knots, and the trawling time varied from 30 to 50 min at each station. The final catch per unit effort (CPUE) was converted into catch  ( weight or quantity ) / net × Hours . An SD204 conductivity temperature depth (CTD) profiler and a handheld GPS locator were used to record the water depth, salinity, latitude, and longitude data of each station ( Table 1 ). After trawling at each station, species identification of teleost fish was conducted based on morphological characteristics, and body weight and body length were measured. The species found in large numbers were sampled for measuring, and those with a small number were all measured.

The water sample collection stations consisted of all 15 trawl survey stations and the other 7 stations ( Figure 1 and Table 1 ). A 2L organic glass hydrophore was used to collect intermediate water, and the water collection at each station was repeated three times. After each collection, the hydrophore was thoroughly washed twice using deionized water. The water samples were filtered through a 0.45μm Glass fiber filter membrane. The filter membrane with enriched eDNA was folded, wrapped with sterile tin foil, stored in the RNA/DNA preservation solution (Vazyme) at -20°C, and transported back to the laboratory. The eDNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen).

The universal primers of marine teleost mitochondrial 12S rRNA gene partial sequence were used for PCR amplification of eDNA (Miya et al., 2015). Primer sequences are as follows: MiFish-F, 5’-GTYGGTAAAWCTCGTGCCAGC and MiFish-R, 5’-CATAGTGGGGTATCTAATCCYAGTTTG (Miya et al., 2015). PCR was repeated three times using the TransStart^® FastPfu DNA Polymerase kit (TransGen). In order to ensure the accuracy and reliability of subsequent data analysis, two conditions should have been met in PCR: 1) Low cycle number amplification should be used as far as possible; 2) The number of cycles amplified for each sample should be consistent. The random selection of representative samples for pre-experiment would ensure that the majority of the samples in the lowest number of cycles could be amplified with an appropriate concentration of products. In this study, the specific PCR conditions were as follows: an initial denaturing at 95°C for 5 min, followed by 30 cycles including 30 s denaturing at 95°C, 30 s annealing at 55°C, 20 s extension at 72°C, and then a final extension at 72°C for 5 min. Triplicate PCR products were pooled, purified, and used for library construction according to Illumina’s genomic DNA library preparation procedure and sequenced using the Illumina MiSeq PE 250 platform (Shanghai Biozeron). All raw sequences were analyzed using QIIME (Caporaso et al., 2010). After filtering out the high-quality reads, the remaining sequences were clustered into operational taxonomic units (OTUs) with a 97% similarity identity cutoff. Classification of the representative sequences of each OTU was performed against the NCBI and the Mitofish database (http://mitofish.aori.u-tokyo.ac.jp/). The corresponding species were then compared with those native to the Bohai Sea and annotated as known native fish as possible.

Due to sample preservation and eDNA extraction failure, among 22 water survey stations, eDNA was successfully extracted from 16 stations, and high-throughput analyses were performed ( Table 2 ). After quality control and screening, the number of sequences obtained from each station varied from 555,030 to 695,056, with an average of 619,009. The average length of the sequences was 167.96 bp ( Table 2 ).

A total of 31 species of marine teleost were detected in 15 survey stations using the trawl survey, and 33 species (including 1 unidentified species, Sebastes spp.) were detected in 16 survey stations using eDNA technology. The number of species detected by both the trawl survey and eDNA was 21 ( Table 3 ). The 21 species included benthic fishes: Amblychaeturichthys hexanema, Cynoglossus joyneri; swimming predatory fishes: Planiliza haematocheilus, Lateolabrax japonicusg; and low economic value fishes: Clupanodon thrissa, S. taty, etc. Combined with the trawl survey and eDNA analysis results, the number of marine teleost species around the Yellow River estuary in June should be no less than 43. Excluding Sebastes spp., eDNA detected 67.74% of trawling species.

Among all 22 survey stations, 10 stations produced both trawling survey data and eDNA high-throughput sequencing results. In this study, the average biomass of each species of these 10 stations was compared with the corresponding average eDNA abundance. The number of species involved in the comparison was 21, i.e., those detected by both the trawl survey and eDNA. The species with the largest catch per hour in the trawl survey was Chaeturichthys stigmatias. Furthermore, this species also had the greatest relative sequence abundance in the water samples ( Table 4 ).

The Pearson’s correlation coefficients between eDNA relative abundance and CPUE (weight and quantity) were both very high (0.993 and 0.993, Table 5 ), and the corresponding determination coefficients (R²) were 0.986 and 0.987, respectively. Furthermore, the Spearman’s correlation coefficients were 0.659 and 0.595, respectively ( Table 5 ). All of the aforementioned correlations were highly significant (p<0.01). There was a clear correlation between biomass and eDNA relative abundance, while there was no significant relationship between relative abundance and average individual body weight ( Table 5 ).

The species composition of 16 survey stations was compared based on the eDNA high-throughput sequencing data of water samples. The results showed that the E2 station had the biggest species composition difference when compared with the other stations, followed by the F1 station, and its β diversity index ranged from 0.326 to 0.555. The β diversity index of the other 15 stations varied from 0.276 to 0.476 ( Figure 2 ). The water depths of E2 and F1 were 4.5m and 3.5m respectively, meaning they were the shallowest of all the 22 stations. Correspondingly, the salinity of the two stations was 25.14 and 15.07, respectively, meaning they also had the lowest salinity levels.

As mentioned above, traditional fishery investigation methods, such as trawling, can obtain first-hand data (species, abundance, and distribution) on fishery species. Although many limitations have emerged in the process of practical application, the trawl survey is the most popular method in fishery investigation at present (Li et al., 2012; Stoeckle et al., 2020; Hongo et al., 2021). Furthermore, due to its inherent uncertainty (Jerde, 2021), eDNA cannot yet replace traditional methods completely; a more accurate fishery resource assessment could be achieved based on a combination of trawl survey and eDNA methods.

In our research, 11 teleost fish species and one genus (Sebastes spp.) were monitored using eDNA but with the trawling method. Conversely, ten species were collected using trawling but with eDNA ( Table 3 ). The discrepancies in fish species assessment between the trawling and eDNA methods should be attributable to the following reasons: the target species is relatively sparse and escaped from trawling, or the special habitats (such as the reef area) where the target species live impeded trawling conduct. In such cases, the target species could be monitored only using eDNA but with the trawling method. For example, the genus Sebastes was mainly distributed in the reef area, where trawling investigation could not be carried out but its eDNA could be monitored. Other potential reasons are that eDNA could have drifted horizontally or vertically with the ocean current (Qian et al., 2021; Wang et al., 2021), or 12S rRNA sequences of some target species may have still been lacking in the databases (such as the NCBI). In these cases, the trawling method could collect living samples but the eDNA result was negative.

Certain measures could be considered to improve the efficiency of eDNA monitoring: 1) Supplementing mitochondrion sequence information on target species in the databases or constructing a special database of native species, thereby increasing the chances of more species being monitored using eDNA technology. 2) In some teleost species within the same genus, such as genus Sebastes in the present study, or genus Takifugu, 12S rRNA sequences are very conservative. In this case, to achieve more accurate discrimination power from eDNA, other genes from the mitochondrion, for example, COI or 18S, or from nuclear genes, might be alternative methods. 3) Due to the “inherent uncertainty” of eDNA (Jerde, 2021), many factors, including temperature, ocean current, water flow, life history, and metabolism, should be considered comprehensively to improve the accuracy of eDNA technology.

Effective fishery management depends on accurate fishery animal biomass assessment. In recent decades, many studies have revealed the correlation between eDNA concentration and biomass or catch, whether using the species-specific eDNA PCR technique (Ct value) or the eDNA high-throughput sequencing (reads abundance) platform. In most cases, the higher the biomass, the higher the eDNA concentration detected (Takahara et al., 2012; Knudsen et al., 2018; Levi et al., 2019; Stoeckle et al., 2020; Wang et al., 2021). The same phenomenon also appeared in the present study; the most abundant OTU was from C. stigmatias, a type of benthic marine fish in the Bohai Sea. Correspondingly, the two measures of biomass (weight and quantity) of C. stigmatias were also ranked the top among all 33 types of fish detected by trawling ( Table 4 ). The second abundant OTU was from A. hexanema, another type of benthic marine fish in the Bohai Sea. Its CPUE (weight) and CPUE (quantity) were both ranked third ( Table 4 ). The CPUE (weight) and CPUE (quantity) of C. joyneri were ranked second and fourth, respectively, and its OTU abundance was ranked fifth ( Table 4 ). In one study, when species-specific eDNA qPCR was involved, significant positive correlations between biomass and eDNA concentration and between CPUE (weight and quantity) and eDNA concentration were detected in Atlantic cod (Gadus morhua) (Salter et al., 2019). In another study in the sea area of New Jersey (USA), in which the eDNA abundance of different marine fishes was calculated using a high-throughput sequencing platform, high positive correlations between OTU abundance and catch were confirmed in more than 70% of the investigated fish species, and these patterns were confirmed in several different fishing seasons (Stoeckle et al., 2020).

In the present study, positive correlations between eDNA abundance and CPUE were confirmed in most of the 21 fish species that were detected by both the trawling method and eDNA ( Table 4 ). Disagreement between eDNA and trawling was found in a few fish species, such as P. haematocheilus, E. japonicas, and Jaydia lineata ( Table 4 ). Interestingly, P. haematocheilus showed a high abundance of eDNA concentration, but few samples were collected ( Table 4 ). A similar result was also revealed in another report (Thomsen et al., 2016), indicating that the trawling method might not be suitable for some large or fast-swimming fish species. Furthermore, although biomass may affect eDNA concentration directly, many other factors may affect it indirectly as well. It is possible that two individuals of the same species and weight can release different eDNA if they feed differently (Klymus et al., 2015). In the future, many factors (such as temperature, development phase, degradation, feeding, and nutrition level) should be considered in order to better understand the correlation between eDNA and biomass.

Using a high-throughput sequencing platform, not only were more fish species able to be monitored than when using the trawling method but the relative abundance or biomass could also be assessed. Furthermore, differences in species composition between the different stations could also be analyzed (Stoeckle et al., 2020). In this study, we also compared the species composition differences among 16 investigation stations based on the composition of marine teleost fish species from the OTU annotation results in each station, from which, stations E2 and F1 displayed the most significant species composition differences when compared with the other 14 stations ( Figure 2 ). Furthermore, we explored the potential environmental factors that might cause these differences: the depths of stations E2 and F1 were the shallowest (4.5 m and 3.5 m, respectively) of the 22 stations, and the lowest salinity levels were also detected in these two stations (25.14‰ and 25.07‰, respectively, Table 1 ). Station E2 lies in the Yellow River estuary and the direction of freshwater flow was to the Southwest and South and into the Bohai Sea. The abundance of freshwater caused sediment deposition and diluted the seawater around stations E2, F1, and F2 (the eDNA sample from station F2 was lost during the experiment) ( Figure 1 ). It is clear that depth and salinity played key roles in affecting the distribution of marine teleost fish distribution in this study. Our findings demonstrate the possibility of using the eDNA technique to trace the effect of circumstances and physical and chemical factors on fish distribution.

It is undeniable, at least at present, that compared with living sample collection, eDNA, to a certain extent, is still not free of uncertainty. However, with continuous improvement and the gradual development of techniques, eDNA could become a robust and precise technique for use in fishery resource research fields (Jerde, 2021).