During the joint survey in the East China Sea from 14 May to 24 May 2019, a total of 44 stations ( Figure 1 ) were designed, and 172 water samples ( Figures 1 and S1 ) were collected using 24 × 20 L Niskin bottles controlled by the SBE 911 Plus CTD system (Sea-Bird Scientific, Bellevue, WA, USA). According to the actual water depth, three layers were defined: the first layer, which was surface water (0–10 m), the second layer (10–50 m), and the third layer (50 m to the bottom). CTD was used to record total depth, water sample depth, temperature, salinity, oxygen, and turbidity, and other environmental variables (chlorophyll, PO₄ ³⁻, NH₄ ⁺, SiO₃ ²⁻, NO₃ ⁻ and NO₂ ⁻) were determined by the staff of the First Institute of Oceanography Ministry of Natural Resources according to the Chinese national standards for each variable (AQSIQ, 1997). Seawater samples were collected using disposable sterile plastic bags, and then 1 L of seawater was immediately filtered using a vacuum pump through a 47-mm glass fiber filter with a pore size of 0.45 μm (Xingya Purification Materials Company, Shanghai, China). To avoid cross-contaminants between sampling sites, all sampling and filtering devices were soaked in 10% sodium hypochlorite solution for more than 10 min before and after each use, cleaned with tap water, and then washed with distilled water. New gloves were replaced after each sample, and 1 L of distilled water was filtered as the negative control on each day of sampling.

After filtration, each filter was folded individually and wrapped in aluminum foil, frozen in liquid nitrogen at −196°C, taken back to the laboratory, and stored at −80°C for further eDNA extraction.

eDNA was extracted from the filters using DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). Inhibitors were removed using OneStep PCR Inhibitor Removal Kit (Zymo Research, Irvine, CA, USA) to improve the quality of eDNA. The purity and concentrations of eDNA were measured by spectrophotometer (NanoDrop 2000, Thermo Fisher Scientific, Waltham, MA, USA).

The Cytb gene sequences of the black sea bream and its taxonomy-related species, such as Acanthopagrus latus, Acanthopagrus butcheri, Diplodus capensis, Diplodus argenteus, Menticirrhus undulatus, Sparus aurata, Sebastes owstoni, and Sebastes taczanowskii, were used to design black sea bream-specific primers and probe ( Figure 2 ). Specific primers and TaqMan probe of the black sea bream were designed using SeqMan, Primer Express 3.0.1, AlleleID 6 software, and primer design tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast). The specificity of the designed primers and probe was verified in silico (primer design tool in National Center for Biotechnology Information (NCBI)) and in vitro (i.e., performing real-time PCR using genomic DNA of black sea bream and its close relatives as the template and DNA extracted from black sea bream pond water as a template).

Black sea bream DNA was extracted by the standard phenol–chloroform–isoamyl alcohol method and amplified using the designed primers. The amplification conditions were as follows: predenaturation at 94°C for 3 min, denaturing at 94°C for 30 s, annealing at 60°C for 30 s and extended at 72°C for 1 min, and 72°C complex extensions for 10 min, for a total of 35 cycles. A 25 μl PCR system was used: 2× San Taq PCR Mix (with Blue Dye) 12.5 μl (Sangon Biotech, Shanghai, China), forward and reverse primers (0.2 μmol/L) each at 1 μl, template DNA 1 μl, and ddH₂O 9.5 μl.

Plasmid standard was prepared by inserting amplified black sea bream PCR products into a pESI-T vector. Then the ligation product was transformed into Escherichia coli DH5α competent cells and cultured at 37°C for 30 min at 180 r/min. Bacterial solution measuring 100 μl was evenly coated on Luria broth (LB) solid culture medium containing ampicillin (100 μg/ml). Inverted culture was performed at 37°C overnight. After blue-white selection, PCR was performed on the bacterial fluid, and the positive results were sent to Sangon Biotech (Shanghai, China) for sequencing. The sequencing results were analyzed and compared for homology. After 24 h of expanded culture, the plasmids were extracted using a high-purity plasmid extraction kit. The purity and concentration of the extracted plasmids were measured by a spectrophotometer (NanoDrop 2000, Thermo Fisher Scientific), diluted to 10⁸ copies/μl of standard, and stored at −80°C.

To avoid contamination, the PCR setup, including preparation and addition of the standard, was performed in one room and qPCR cycling in another. Real-time PCR was performed with a 7300Plus Real-Time PCR System (ABI, Foster City, CA, USA). The amplification conditions were as follows: denaturation at 94°C for 5 s, and annealing, extension, and data acquisition at 60°C for 30 s, for a total of 40 cycles. The reaction system (20 μl) contained TaqMan™ Fast Advanced Master Mix (Thermo Fisher Scientific) 10 μl, each primer (10 μM) 0.4 μl, TaqMan probe (10 μM) 0.2 μl, 4 μl eDNA extracted, and PCR-grade water 5 μl (Sangon, Shanghai, China). Each experiment was repeated three times, and three negative controls (pure water as the template) were set. Plasmid DNA was diluted from 10⁸ copies/μl in a 10-fold concentration gradient to 10¹ copies/μl as a standard for quantitative PCR. Series dilutions of plasmid standards were used to determine the limit of detection (LOD) and limit of quantification (LOQ) of the whole method.

All statistical analyses were performed using R software version 3.5.3. Since the data between water layers and stations did not meet the normality and homogeneity of variance, statistical differences in black sea bream eDNA concentration among the three layers and stations were assessed using the Kruskal–Wallis test. The average eDNA concentration was calculated for each water layer. It is calculated as follows: the average for each layer is the sum of all eDNA concentrations in that water layer divided by the number of stations.

The relationship between environmental variables and log-transformed black sea bream eDNA concentrations (i.e. Ln(eDNA + 1)) was calculated using the 'pheatmap' package.

Multiple linear regression (MLR) was used to determine the best fitting models for black sea bream eDNA concentrations based on environmental variables. A preliminary analysis was first carried out to ensure that the assumptions of normality, linearity, and homoscedasticity were not violated. The eDNA concentrations were logarithm transformed to conform to a normal distribution prior to analysis. Within the MLR model, the 12 environmental variables (i.e., total depth, water sampling depth, temperature, salinity, oxygen, turbidity, chlorophyll, PO₄ ³⁻, NH₄ ⁺, SiO₃ ²⁻, NO₃ ⁻, and NO₂ ⁻) were included as independent variables, and Ln([eDNA]+1) was included as the dependent variables. The MLR was carried out utilizing all subset regression methods of the ‘Leaps’ package. The adjusted Akaike information criterion (AIC) and the ‘gvlma’ package were used to validate the MLR model. Co-linearity between predictive variables was assessed using variance inflation factors (VIFs), and covariates with VIF greater than 5 were removed.

The sampling stations were visualized with ArcGIS 10.1 software. Gridding of black sea bream eDNA concentrations data at each sampling station was performed using Kriging interpolation and culling of gridded data from the study area using mask extraction.

The spatial distribution of black sea bream was visualized with Ocean Data view 4 software. This was performed as follows: black sea bream log transformed eDNA concentrations at each sampling site were gridded by the DIVA gridding method, adjusting the scale lengths to cover the entire study area. Other parameters were selected as default values and applied to all graphs.

A set of black sea bream-specific primers and a probe were developed to amplify a 146-bp sequence based on Cytb. The sequence was as follows: black sea bream-cytb-F: CTTCTCGGTCTCTGCTTAATCTCC, black sea bream-cytb-R GCGTGGAGGTTTCGGATTAG, black sea bream-cytb-P: ACTTCCGACATTGCCAC ( Figure 2 ). The initial concentration of the plasmid standard was 2.78 × 10⁸. LOD and LOQ were 27.80 and 84.24, respectively. In the results of qPCR, the Ct value of the plasmid standard was 13.3–32. The correlation coefficient, R², of the curve was 99.99%. Amplification efficiency was 107.074%. Within 40 cycles, only samples from black sea bream DNA as the template by using a detection limit of 27.8 were detected, while other species and negative controls were not detected. The results showed a good linear relationship within the standard range of dilution plasmid, and the standard curve established correctly reflected the amplification of Cytb gene in black sea bream.

The black sea bream eDNA was detected in 35 samples (20.35%) at only eight stations (18.18%) of all the sites sampled. There is no significant difference in the black sea bream eDNA concentration among 44 sampling sites, indicating no significant difference in the horizontal distribution of black sea bream eDNA. The average concentration of eDNA was 1.238 ± 1.359 copies/L. Among the three horizontal water layers, the eDNA concentration of Layer 2 (1.481 ± 2.426 copies/L) was higher than that of the two other water layers (Layer 1, 0.883 ± 1.655 copies/L; water Layer 3, 1.257 ± 1.878 copies/L). There was no significant difference between the different water layers. The eDNA hotspots were mainly concentrated in the Yangtze Estuary fishing ground (S00-2), the Dasha fishing ground (S00-3 and S00-4), and the Zhoushan fishing ground (S01-2 and S01-3) ( Figures 3 ).

Although no significant difference was observed in the vertical (difference between layers) black sea bream eDNA concentration distribution either, a higher concentration of eDNA was detected in water depths approximately 30–40 m, such as S00-1, Layer 2; S01-2, Layer 2; S01-3, Layer 2; S02-4, Layer 2; and S03-4, Layer 2 ( Figure 4 ). For different longitudes, the eDNA hotspots were mainly concentrated at 29°N–33°N, within which the average concentration was 19.7541 copies/L. The latitude section shows that the eDNA hotspots were mainly in 122°E–125°E, and a significant difference existed between the sampling points in this range and those outside (p = 0.004212 < 0.01).

The correlation analysis of the 12 environmental variables at the measured sites revealed that black sea bream eDNA presence or absence and concentration are positively correlated with the temperature of sampling stations at a significant level (p < 0.05), as shown in Figure 5 .

In MLR analysis, environmental factors and eDNA concentration were predicted as follows: ln[(eDNA) + 1] =0.0012 + 0.87 × Temperature + 1.15 × NO₃ ⁻+ 1.37 × Salinity + 0.15 × PO₄ ³⁻ (PO₄ ³⁻, p <0.01,Temperature, p < 0.05, NO₃ ⁻, p < 0.05, Salinity, p < 0.05).

This study assessed the potential distribution of black sea bream in the East China Sea by measuring the eDNA presence/absence and concentrations at different sampling stations. From the perspective of horizontal distribution, black sea breams are mainly distributed in the Yangtze Estuary fishing ground, the Dasha fishing ground, and the Zhoushan fishing ground. The largest estuary in China and the western Pacific, the Yangtze estuary, has large seasonal variations in environmental factors such as temperature and salinity and has high primary productivity (Zhang et al., 2017). A total of 1,390 islands and more than 3,850 reefs are distributed in the sea area of the Zhoushan Archipelago. Many sessile organisms inhabit these rocky reefs (Saito et al., 2008; Bissegger et al., 2014). The excellent natural conditions of these sea areas provide an ideal habitat for black sea bream. Therefore, black sea bream is concentrated in the areas of the Yangtze River estuary and the Zhoushan Archipelago (Yamada et al., 2007; Liang et al., 2010).

Monitoring eDNA concentration enables predicting a target species’ potential distribution in the target area, reducing workload and research costs, and enabling more effective ecological monitoring (Erickson et al., 2016; Buxton et al., 2017). The spawning seasons of black sea bream populations vary between different sea areas. Black sea bream spawning times are as follows: from December to mid-March of the following year in Guangdong and Taiwan areas, from November to March of the following year in Hong Kong and its adjacent area (Law and Sadovy de Mitcheson, 2017), mid-March to early May in Fujian (Urbatzka et al., 2012), from early May to early June in Shandong (Zhang et al., 1980), and from April to May in Jiangsu (Cai, 1983). The water samples were collected in mid-May, close to the spawning time of black sea bream in Jiangsu, and the results showed that the black sea bream eDNA signals were high in the surrounding area of Jiangsu, which strongly suggests the presence of spawning fish. Previous studies have confirmed the ability of aquatic species spawning ground detection using eDNA (Takeuchi et al., 2019; Wang et al., 2020), for example, the use of environmental DNA monitoring to confirm the results of changes in spawning grounds of small yellow croaker (Larimichthys polyactis) that were first identified by traditional methods (Wang et al., 2020). The mystery spawning region of the Japanese eel Anguilla japonica was eventually discovered and validated by the eDNA approach (Takeuchi et al., 2019). In this study, environmental DNA hotspots coincided in time and extent, with black sea bream spawning grounds documented by traditional surveys (Yamada et al., 2007). We therefore speculated that the areas around S00-2 and S00-3 were the spawning sites of black sea bream in the sea near Jiangsu. Although this speculation will require more detailed field recruitment work, the use of the eDNA approach has greatly improved the efficiency of species monitoring.

The vertical distribution of eDNA depends on the location of organisms and their behavior (Murakami et al., 2019b). By investigating and researching the vertical distribution, species distribution and habitat preferences can be more comprehensively mastered. For example, changes in the eDNA concentration in the surface and deep layers during dawn and dusk can be used to estimate the times when organisms begin and end daily migrations (Govindarajan et al., 2021). eDNA concentration changes reflect the fish seasonal variation in thermal stratification and habitat preferences (Littlefair et al., 2021). Studies have shown that even when organisms are migrating, there is a residual eDNA signal at the major inhabited by the species and that this signal does not disappear completely (Allan et al., 2021). The black sea bream eDNA is mainly concentrated in the shallow zone of 20–40 m, which may be related to the fact that black sea bream live and breed mainly in this depth range (Wang et al., 2015). Black sea bream eggs are pelagic, and the planktonic life of the eggs and larvae estimated at about 20 days suggests that black sea bream can spread far from the current (Jean et al., 1998). In the sea area around the Yangtze Estuary, the distribution pattern of early fish larvae is more in the north and less in the south (Li et al., 2020). The high-density area of larvae is often on the offshore side of the spawning ground (Jiang et al., 2020), which can explain the high eDNA concentrations in surface water detected near the estuary of the Yangtze River. More eDNA was detected in the surface layer from the Mindong fishing ground (S05-5 and S05-6), which may be under the influence of the Taiwan warm current from March to May (Lu and Lee, 2014). Although the results of eDNA monitoring can be explained by the habits of black sea bream, the lack of traditional resource survey results prevents simultaneous validation of the correlation between environmental DNA and black sea bream resources. In addition, eDNA is unable to distinguish between the size, age, and survival status of target organisms (Nevers et al., 2018; Deutschmann et al., 2019). eDNA can also be dispersed and degraded by complex environmental factors in the ocean, such as currents (Yamamoto et al., 2016; Murakami et al., 2019a), water temperature (Jane et al., 2015; Lacoursière-Roussel et al., 2016b), pH, and UV (Hansen et al., 2018). However, research on these issues has focused more on closed and semi-closed environments. Fewer studies have been conducted on large open areas (Yamamoto et al., 2016; Murakami et al., 2019a). Therefore, given the obvious knowledge gaps that still exist regarding the longevity and influencing factors of eDNA signaling, eDNA studies should be seen as complementary rather than a ‘stand-alone tool’, especially in the flowing sea area. More studies are needed to further improve our knowledge of the parameters and processes that could potentially affect the eDNA distribution in the ocean.

Among the three water layers, the highest average concentration of eDNA was found in the second layer. This difference is inextricably related to halocline, pycnocline, and eDNA degradation (Hansen et al., 2018; Takasu et al., 2019; Jeunen et al., 2020). In addition, water temperature and bait distribution are influential factors affecting the vertical distribution of fish (Gibson, 1994). The results of water sample display at a sampling site may not necessarily represent the biological condition of the site (Nakagawa et al., 2018). Considering target species habitat preference and the water column environment differences, both the horizontal and vertical distributions of organisms need to be taken into account in future surveys of other marine organisms to better assess fishery resources.

The survival and growth of fish are inextricably linked to the marine environment. Thus, environmental factors play pivotal factors in the investigation of fishery resources (Mukherjee et al., 2012). In the predicted regression equation, in addition to the previously mentioned temperature factors that would influence black sea bream eDNA concentrations, PO₄ ³⁻, NO₃ ⁻, and salinity may influence eDNA concentration too. The concentration and the presence or absence of eDNA of black sea bream are also significantly correlated with temperature (p < 0.05). The black sea bream’s most suitable growing temperature is 17°C–25°C, and the spawning temperature is 16°C–21°C (Jiang et al., 2012). The average temperature, 15.51°C, at the sampling stations near Jiangsu is close to the spawning temperature, resulting in more eDNA being detected. The concentrations of PO₄ ³⁻, in a few ranges where eDNA hotspots showed up, are between 0 and 0.5 mg/L, while the concentrations of NO₃ ⁻ are within the extent of 0–12 mg/L. The concentrations of both PO₄ ³⁻ and NO₃ ⁻ are at normal levels. A particular concentration of PO₄ ³⁻ and NO₃ ⁻ is conducive to algae growth, increasing primary productivity and providing bait directly or indirectly (Sarre et al., 2000). Salinity is also an essential factor affecting the distribution and abundance of black sea bream. Although this species is a broad-salinity fish, only the sea area with high salinity can promote its gonad maturity and spawning (Murata et al., 1997). The spatial distribution of fish stocks and their relationship with the environment has always been a research hot spot in fishery ecology (Hu et al., 2020). Consequently, we propose that in future surveys of fishery resources in the open ocean, consideration should be given to inferring the distribution patterns under different environmental conditions by predicting the regression equations between environmental factors and fish, so as to obtain the spatial distribution characteristics and influencing factors of target fish and conduct surveys with more specificity. This proposal will provide theoretical reference and support for the development of management strategies, such as the sampling design of fish communities and the establishment and optimization of fish nature reserves.

In the present context, the data obtained from quantitative eDNA monitoring methods remain questionable due to several complex environmental factors (Hinlo et al., 2018; Tillotson et al., 2018; Antognazza et al., 2021). Thus, despite the relatively high eDNA detection rates in this study, information on the spatial resolution of these detections (distribution and biomass) is inconclusive. More studies are needed in the future to validate and complement the data provided by eDNA methods, such as traditional fishery surveys for confirmation at eDNA signal hotspots.