Species Composition and Assemblages of Ichthyoplankton in Sansha Bay, Fujian Province, China

Sansha Bay (26.40−27.00°N, 119.50−120.20°E) is a typical semi-enclosed bay, located in northern Fujian Province, China, and adjacent to the East China Sea. The ichthyoplankton species composition and assemblage structure were investigated based on monthly sampling at 25 stations in April−September 2019, covering the important spring and summer spawning seasons in the region. Sampling was conducted in the first 3−5 days of the full moon or new moon phases using a standard plankton net through horizontal and vertical tows during daytime. In total, 25,819 ichthyoplankton samples were collected, of which 25,449 samples (i.e., 24,757 eggs and 692 larvae) were from horizontal tows. For horizontal tow samples, the ichthyoplankton were classified into 58 taxa in 15 orders and 23 families with a combination of external morphology and DNA barcoding analyses, from pelagic to demersal and benthic species. The dominant order was the Gobiiformes, including 23 species (39.7% of all species). The dominant taxa, in terms of relative abundance and frequency of occurrence, consisted of commercially important fishes, such as Setipinna tenuifilis (Valenciennes, 1848) (Engraulidae), Epinephelus akaara (Temminck and Schlegel, 1842) (Serraenidae), Collichthys lucidus (Richardson, 1844), Nibea albiflora (Richardson, 1846) (Sciaenidae), Acanthopagrus schlegelii (Bleeker, 1854), and Pagrus major (Temminck and Schlegel, 1843) (Sparidae), accounting for 78.9% of the horizontal tow samples. Low-valued and small-sized fishes, such as Stolephorus commersonnii Lacepède, 1803 (Engraulidae), Solea ovata Richardson, 1846 (Soleidae), Nuchequula nuchalis (Temminck and Schlegel, 1845), and Photopectoralis bindus (Valenciennes, 1835) (Leiognathidae), were also dominant species, accounting for 11.4% of the horizontal tow samples. The ichthyoplankton assemblage was categorized into five different temporal assemblages based on the cluster and nonmetric multidimensional scaling analysis, namely, April, May, June, July, and August−September (ANOSIM, Global R = 0.656, p < 0.01) with the highest density and richness of ichthyoplankton occurred in May. The spatial distribution pattern showed that the high density (ind./m3) of ichthyoplankton occurred mainly in S12–S25 in Guanjingyang and along the Dongchong Peninsula coastline into Dongwuyang, while low density occurred mainly in S01–S11 in the northwest waters of Sandu Island (ANOVA, F = 8.270, p < 0.05). Temperature, salinity, and chlorophyll a were key factors structuring the ichthyoplankton assemblages in Sansha Bay. In addition, this study revealed the changes of the ichthyoplankton composition, density, and spatial distribution in Sansha Bay over the past three decades.

Sansha Bay (26.40−27.00 • N, 119.50−120.20 • E) is a typical semi-enclosed bay, located in northern Fujian Province, China, and adjacent to the East China Sea. The ichthyoplankton species composition and assemblage structure were investigated based on monthly sampling at 25 stations in April−September 2019, covering the important spring and summer spawning seasons in the region. Sampling was conducted in the first 3−5 days of the full moon or new moon phases using a standard plankton net through horizontal and vertical tows during daytime. In total, 25,819 ichthyoplankton samples were collected, of which 25,449 samples (i.e., 24,757 eggs and 692 larvae) were from horizontal tows. For horizontal tow samples, the ichthyoplankton were classified into 58 taxa in 15 orders and 23 families with a combination of external morphology and DNA barcoding analyses, from pelagic to demersal and benthic species. The dominant order was the Gobiiformes, including 23 species (39.7% of all species). The dominant taxa, in terms of relative abundance and frequency of occurrence, consisted of commercially important fishes, such as Setipinna tenuifilis (Valenciennes, 1848) (Engraulidae), Epinephelus akaara (Temminck and Schlegel, 1842) (Serraenidae), Collichthys lucidus (Richardson, 1844), Nibea albiflora (Richardson, 1846) (Sciaenidae), Acanthopagrus schlegelii (Bleeker, 1854), and Pagrus major (Temminck and Schlegel, 1843) (Sparidae), accounting for 78.9% of the horizontal tow samples. Low-valued and small-sized fishes, such as Stolephorus commersonnii Lacepède, 1803 (Engraulidae), Solea ovata Richardson, 1846 (Soleidae), Nuchequula nuchalis (Temminck and Schlegel, 1845), and Photopectoralis bindus (Valenciennes, 1835) (Leiognathidae), were also dominant species, accounting for 11.4% of the horizontal tow samples. The ichthyoplankton assemblage was categorized into five different temporal assemblages based on the cluster and nonmetric multidimensional scaling analysis, namely, April, May, June, July, and August−September (ANOSIM, Global R = 0.656, p < 0.01) with the highest density and richness of ichthyoplankton occurred in May. The spatial distribution pattern showed that the high density (ind./m 3 ) of ichthyoplankton occurred mainly in S12-S25 in Guanjingyang and along the Dongchong Peninsula coastline into Dongwuyang, while low density occurred mainly in S01-S11 in

INTRODUCTION
Fishes have various life stages, and different developmental stage requires different foods, habitats, and environmental factors (Rijnsdorp et al., 2009;Hamre et al., 2013;Petitgas et al., 2013;Madeira et al., 2020). For example, ripe fishes are sensitive to the hydrological conditions of the spawning grounds, such as current, tide, water temperature, salinity and dissolved oxygen, and prey availability and preference (MacKenzie et al., 1996;Costa et al., 2002). Ichthyoplankton include fish eggs and larvae, belonging to the early life stages. Due to their absence of or weakness on independent swimming capabilities, ichthyoplankton have a drifting nature with currents and tides (Lechner et al., 2016;Downie et al., 2020). Understanding the distribution patterns of different life stages and exploring the abiotic and biotic factors influencing the distribution patterns are crucial for assessing fish recruitment and stock restoration (Costa et al., 2002;Santos et al., 2017). Achieving the sustainability of fishery resources depends highly on the abundance and survival of fish eggs and larvae (Oeberst et al., 2009;Llopiz et al., 2014).
The identification of unambiguous species from fish eggs and larvae is challenging. It is mainly because of the high diversity of species of fishes and the limitation of morphological observation on different stages of eggs and larvae before metamorphosis (Li et al., 2014;Zhang et al., 2015;Hsieh et al., 2016). For example, eggs and larvae of Sciaenidae are difficult to identify to species level through morphological features (Zhang et al., 1985;Zhan et al., 2016), besides the challenges for identification of adult sciaenids (Chu et al., 1963). DNA barcode, which delimits species using a molecular marker (e.g., mitochondrial cytochrome c oxidase subunit I gene, COI), has been widely applied to facilitate the identification in all stages from egg to adult of fishes (Ward et al., 2005;Ko et al., 2013;Becker et al., 2015;Harada et al., 2015;Hou et al., 2021a). In recent studies, Kerr et al. (2020) successfully identified 564 fish eggs to 89 taxa collected from northwestern Cuba and across the Florida Straits using DNA barcode. Hou et al. (2021a) successfully identified 931 fish eggs and 229 larvae to 75 taxa collected from the Pearl River Estuary of China using a DNA barcode. However, identifying ichthyoplankton correctly through DNA barcoding relies highly on a powerful and accurate database. It is necessary, but still challenging, to build up a DNA barcoding library in a specific area.
The transitional region between the East China Sea and the South China Sea through the Taiwan Strait, Fujian Province, has a highly sinuated coastline forming many bays with some typically semi-closed by peninsulas (Yu et al., 1988). In terms of volumes of the marine capture fisheries, Fujian Province ranks third in China, contributing to approximately 15.6% of the national total marine capture volumes or 16.8% of the national total marine fish capture volumes in the past decade, 2010−2019 (MOA, 2010(MOA, -2018MARA, 2019). Sansha Bay (26.40−27.00 • N, 119.50−120.20 • E) is located in Fujian Province ( Figure 1A) and connects with the East China Sea through the narrow Dongchong Channel of about 3 km in width ( Figure 1B). This typical semi-enclosed bay is well-known because it was one of the 12 traditional spawning grounds along Chinese coastal waters for the large yellow croaker Larimichthys crocea (Richardson, 1846) (Sciaenidae) before the 1990s, a commercially important species that formed large spawning aggregations (Chu, 1985;Liu and Sadovy de Mitcheson, 2008). With the mariculture and industrial development since the 1990s, marine pollution has been exacerbated in Sansha Bay (Wang et al., 2019). The current status of Sansha Bay as L. crocea spawning ground function is unclear.
In Sansha Bay, species composition and assemblages of ichthyoplankton, which are essential for the estimation of spawning sites, reproductive seasons, and conservative strategies (Parrish et al., 1981;Bailey and Houde, 1989;Oliveira and Ferreira, 2008), remain insufficient attention. Species composition and temporal and spatial distribution of fish eggs and larvae were investigated only four times (i.e., 1990, 2007, 2008, and 2010) historically, mainly seasonally, and the samples were identified by external morphology only (Dai, 2006;Wang et al., 2010;Shen, 2011;Xu, 2018), which may lead to the missing of certain species and underestimate the diversity of the real species. The results of these previous studies showed that Sansha Bay was important for fish reproduction, and fish egg and larva density (ind./m 3 ) were high in spring and summer. Moreover, the ichthyoplankton abundance of traditionally important fishery species showed dramatic declines over years. It merits further monthly and accurate studies with multiple identification methods and a high spatial coverage to evaluate the current status and long-term dynamics of ichthyoplankton assemblages in Sansha Bay.
In this study, we designed 25 stations throughout Sansha Bay for the ichthyoplankton collection. We conducted sampling monthly from April to September in 2019 covering the peak spawning seasons (i.e., spring and summer) mentioned earlier and applied DNA barcoding technique for species identification after external morphology examination, aiming to obtain a relative accurate species composition, to provide detailed characteristics of the ichthyoplankton assemblages, and to analyze the influence of abiotic and biotic variables on the distribution of ichthyoplankton in Sansha Bay. The results can be used to evaluate the current status of this traditional spawning ground and to compare the changes of Frontiers in Marine Science | www.frontiersin.org ichthyoplankton composition and spatial distribution over the past three decades in Sansha Bay.

Ichthyoplankton Collection
Ichthyoplankton samples were collected from 25 stations (S01−S25) (26.50−26.85 • N, 119.52−120.09 • E) in Sansha Bay, Fujian Province of China ( Figure 1B). The sampling was conducted monthly from April to September 2019 within the first 3−5 days of the full moon or new moon phases during daytime.
The standard plankton net used for horizontal tows is 80 cm in net opening diameter (i.e., 0.5 m 2 net opening area) and 505 µm in net mesh size (GB/T 12763.6, 2007). At each station, the net was towed horizontally right below the sea surface with the speed of approximately 2 km for 5−10 min.
The standard plankton net used for vertical tows is 50 cm in net opening diameter (i.e., 0.2 m 2 net opening area) and 505 µm in net mesh size (GB/T 12763.6, 2007).
On board, ichthyoplankton samples collected were fixed with the formaldehyde concentration of 5%. The samples were further transferred to 95% ethanol solution within 12 h for the purpose of DNA barcoding.

Collection of Abiotic and Biotic Factors
Sea surface temperature, salinity, dissolved oxygen, and pH were measured using YSI EXO2 Multiparameter Sonde during plankton sampling at each station.
Chlorophyll a (Chla) samples were collected and analyzed in the laboratory through fluorospectrophotometry within 24 h. For a sampling of vertical phytoplankton, the standard plankton net used is 37 cm in net opening diameter (i.e., 0.1 m 2 net opening area) and 77 µm in net mesh size (GB/T 12763.6, 2007). On board, phytoplankton samples were fixed in Lugol's iodine solution with a concentration of 10%. For a sampling of vertical zooplankton, the standard plankton net used is the same as the vertical ichthyoplankton collection above. On board, zooplankton samples were fixed in the formaldehyde solution with a concentration of 5%.

Ichthyoplankton and Plankton Identification
In the laboratory, all fish eggs and larvae collected were sorted and pooled under a dissecting microscope (Zeiss Stemi 2000-C, Germany) according to the external morphological features (Zhang et al., 1985;Shao et al., 2001;Wan and Zhang, 2016). For each egg or larva with a specific morphological feature, photo images were collected using the Leica M165FC Fluorescent Stereo Microscope and Canon EOS 600D Digital SLR Camera, and one or more samples were selected for DNA extraction and barcode COI fragment obtaining. Sizes of eggs and larvae were not measured because they were fixed and dehydrated in 95% ethanol solution for at least 2 months before sorting.
Genomic DNA was extracted using the Takara MiniBEST Universal Genomic DNA Extraction Kit Ver. 5.0 according to manufacturer specifications and further used with no dilution for amplification and sequencing. A PCR was conducted, and the COI fragment was amplified using three universal COI gene primer pairs (i.e., LCO1490−HCO2198, FishF1−FishR1, and FishF2−FishR2) (Folmer et al., 1994;Ward et al., 2005) to obtain the target COI gene fragment (about 630 bp). For each PCR product, sequencing was performed using both corresponding forward and reverse primers that are used for PCR (Sangon Biotech, China).
Raw sequences were assembled using Bioedit 7.2. Assembled sequences of fish eggs and larvae were identified using reference sequences published in GenBank (NCBI) 1 or by the Barcode of Life Data System. 2 To assign a sequence to a species, genus, or family level, we required that the sequence matches with the database by at least 99.60, 91.29, or 85.48%, respectively (Hou et al., 2017). All sequences were further confirmed by ML and NJ phylogenetic trees with sequences from the GenBank and BOLD system using MEGA 6.06 with 1,000 bootstrap pseudoreplications. Categorization of species taxonomy was followed as suggested by Nelson et al. (2016).

Plankton Calculations
In the laboratory, the densities of ichthyoplankton (i.e., eggs and larvae) from horizontal tows were standardized: where A is the density of ichthyoplankton (ind./m 3 ), N is the number of ichthyoplankton (individuals) per plankton net, and V is the seawater volumes filtered into the plankton net (m 3 ).
V was calculated by a flow meter (HYDRO-BIOS) attached to the center of the plankton net: where R is the number of revolution noted from HYDRO-BIOS, 0.3R means that when the flow meter rotates one number, the net moves 0.3 m based on HYDRO-BIOS product manual, and S is the net opening area (m 2 ).
The dominant species of ichthyoplankton (i.e., eggs and larvae together) from horizontal tows were determined using the Index of Relative Importance (IRI) (Pinkas et al., 1971;Zhu et al., 2002): where N% and F% are relative abundance and frequency of occurrence of a specific species, respectively. Species with IRI ≥ 0.02 were considered to be a dominant species.
The three diversity indices, namely, the Shannon-Wiener diversity index (H ), Pielou's evenness index (J ), and Margalef 's richness index (d), were calculated (Margalef, 1958;Shannon and Wiener, 1963;Pielou, 1966) to evaluate the current status of ichthyoplankton in Sansha Bay: where S is the total number of ichthyoplankton species, N is the total number of individuals, and P i is the proportion of the number of individuals for a specific ichthyoplankton species to the total number of individuals.
In addition, densities of phytoplankton and zooplankton from vertical tows were also calculated as for ichthyoplankton aforementioned: where A is the density of phytoplankton (cells/m 3 ) or zooplankton (ind./m 3 ), N is the number of phytoplankton (cells) or zooplankton (individuals) per plankton net, and V is the seawater volumes filtered into the plankton net (m 3 ).

Data Analysis
We used only data from horizontal tows for further analyses for at least two reasons. First, the number of eggs and larvae and diversity of ichthyoplankton species from horizontal tows were much higher than those of vertical tows in sampling during every month. Second, the depths of sampling stations varied largely from as shallow as 1 m to about 40 m. Therefore, the seawater volumes filtered into the vertical tow nets varied largely that it subsequently influenced the number of eggs and larvae and the diversity of ichthyoplankton species collected. The vertical tow samples were only used for species diversity supplements. For the data analysis below and specific calculation above (e.g., IRI), we combined the number of egg and larva samples together Huang et al., 2017;Hou et al., 2021a). A short hatchery period (i.e., hours to several days) and limited independent swimming capability of larvae were also the considerations for the combination of egg and larva data for analyses (Woynarovich and Horváth, 1984;Pauly and Pullin, 1988;Lechner et al., 2016;Downie et al., 2020).
The one-way ANOVA test was conducted with the IBM SPSS Statistics version 16.0 for revealing the significance of the monthly changes of environmental factors, the diversity indices, the spatial variation of average ichthyoplankton density, and the variation of average species richness among 25 stations; when necessary (ANOVA, p < 0.05, n > 2), a post hoc least significant difference (LSD) test was added. The cluster analysis and nonmetric multidimensional scaling (nMDS) based on the Bray-Curtis similarity and square root average density of ichthyoplankton were performed to clarify the temporal and spatial distribution patterns of the assemblages. After the cluster analysis, the tests for identifying the discrimination between two observed sample clusters (SIMPEROF) and the similarities (ANOSIM) were performed using PRIMER version 6.0 to assess the significance and discrimination between sample clusters (Clarke and Gorley, 2006).
The detrended correspondence analysis (DCA) and canonical correlation analysis (CCA) or redundancy analysis (RDA) were performed using CANOCO version 5.0 to analyze the constrained relationships between environmental factors and ichthyoplankton communities. Before the analysis, DCA was performed to determine a suitable response according to the maximum gradient of length (>4, CCA; 3−4, CCA or RDA; <3, RDA) (Šmilauer and Lepš, 2014). Because the maximum gradient length of DCA was 6.34 (>4), CCA was subsequently selected in this study. Only species that occurred in >10% of the catches, based on all species density, were included in the analysis to reduce the effect of rare species. Data of species density were log (10,000x + 1) transformed to minimize the dominant effect of some species (Clarke et al., 2014). A Monte Carlo permutation test with 999 permutations was performed to confirm the key factors significantly affecting the assemblages (p ≤ 0.05; Arora and Mehra, 2009;Hou et al., 2021a). The variance inflation factor (VIF) was conducted to examine the collinearity between independent environment variables. VIF < 10 indicated that the environmental variables were not collinear and were fit for CCA (Graham, 2003;Huang et al., 2017). The importance of the environmental variables was measured by interset correlation coefficients, when the value of interset correlation coefficients ≥| ±0.4| variables were regarded conservatively as biologically important (Rakocinski et al., 1996;Ramos et al., 2017).

Ichthyoplankton Composition
A total of 25,449 ichthyoplankton samples, including 24,757 fish eggs and 692 larvae, were collected from horizontal tows. A total of 370 ichthyoplankton samples, including 315 eggs and 55 larvae, were collected from vertical tows. The numbers of eggs and larvae from horizontal tows were much higher than those of vertical tows.
According to different morphological features, 619 fish eggs and 192 larvae were selected for photography and further DNA extraction, from which high-quality COI fragments were obtained from 335 fish eggs (54.1%) and 117 larvae (60.9%). In total, 60 taxa were identified, 52 at species level, including 26 from eggs and 43 from larvae (Supplementary Table 1).  Figure 1). Additionally, four larvae were identified at the genus level (i.e., Pseudogobius and Rhinogobius in Gobiiformes, Hyporhamphus in Beloniformes, and Hippichthys in Syngnathiformes), and three larvae at the family level (i.e., Gobiidae in Gobiiformes and Blenniidae in Blenniformes).
From vertical tows, 31 taxa were identified from 9 orders, 13 families with 27 at the species level, 3 at the genus level, and 1 at the family level (Supplementary Table 2). Among these, 29 taxa were also recorded in horizontal tow samples, and only two species [i.e., Acentrogobius sp. and Lateolabrax japonicus (Cuvier, 1828)] were only found in vertical tow samples. Eggs and larvae of two alien species, the red drum Sciaenops ocellatus (Linnaeus, 1766) (Sciaenidae) in Acanthuriformes and the gilthead seabream Sparus aurata (Linnaeus, 1758) (Sparidae) in Spariformes, were collected in horizontal tows, S. ocellatus eggs in August and September and S. aurata larvae in June (Supplementary Table 1). Eggs and larvae of L. crocea, a critically endangered species in the IUCN Red List, were collected in April−June in horizontal and vertical tows (Supplementary Tables 1, 2).
Ichthyoplankton from Engraulidae (two species), Soleidae (one species), Serranidae (one species), Leiognathidae (two species), Sciaenidae (two species), and Sparidae (two species) were dominant (IRI ≥ 0.02; The three diversity indices showed significant differences from April to September (H : F = 10.39, p < 0.05; J : F = 2.54, p < 0.05; d: F = 16.14, p < 0.05; Figure 2). The highest H and d were found in May, while the lowest H and d were found in September; J in August was significantly higher than the other months.
The ichthyoplankton were classified into five representative assemblages in the Bray-Curtis similarity at 36.79% (p = 0.001), namely, April, May, June, July, and August−September ( Figure 4A). The results of nMDS ( Figure 4B) are consistent with those of the cluster analysis at the stress value < 0.2,   Table 3). The ichthyoplankton assemblages of April, May, June, and July differed significantly from August and September, and no significant difference was found between August and September.
Although the cluster analysis of spatial assemblage was not significant (SIMPEROF, p > 0.05), ichthyoplankton density (total taxa) ( Figure 5) and richness (Figure 6) varied by station and month. The highest density and species richness occurred in May, from May-S22 and May-S01, respectively. High density occurred mainly in S12-S25 in Guanjingyang and along the Dongchong Peninsula coastline (i.e., southeast Sansha Bay) into Dongwuyang (i.e., northeast Sansha Bay) in April−September, except for August, and low density occurred mainly in S01-S11 in northwest waters of Sandu Island in April−September (ANOVA, F = 8.270, p < 0.05). Although the spatial patterns for the species richness were not clear, it showed that around Sandu Island, especially S02, S03, S05, and S06, had a relatively high richness (ANOVA, F = 7.953, p < 0.05) during most of the sampling period.

Relationship Between Ichthyoplankton Assemblages and Abiotic and Biotic Factors
The variance inflation factor (VIF) test and Monte Carlo permutation test indicated that abiotic factors (i.e., T, S, DO, and pH) and biotic factors (i.e., Chla, Phy, and Zoo) were not collinear (VIF < 10) and contributed significantly to explain ichthyoplankton assemblage structure (p < 0.05; Table 3).
The sum of all canonical eigenvalues occupied 18.12% of the sum of all eigenvalues. The cumulative percentage variance of species was 15.79%, and cumulative percentage variance of species environment was 87.12% ( Table 4). The first two CCA axes (i.e., axes 1 and 2) explained 73.65% of the cumulative percentage variance of species and 73.67% of the cumulative percentage variance of species environment.

Ichthyoplankton Composition and Variation
We reported the first results combining the external morphology and DNA barcode technique to identify ichthyoplankton samples collected from April to September 2019 in Sansha Bay, which made the species richness increase significantly. In total, 60 taxa were identified from horizontal and vertical tows in this study, significantly higher than previous surveys (23−46 species in 1990−2010) from the same sampling area ( Table 5). The high ichthyoplankton richness in this study was mainly from Gobiidae (21 species) and contributed to the application of COI gene fragment analysis with the high detection rate at the species level, rising from 39.1 to 74. 2% (1990−2010) to more than 87.10% (2019) ( Table 5).
The development of local or regional fish species COI barcodes library is highly recommended (Hou et al., 2021b). We demonstrated the importance of the accuracy and integrity of COI barcodes database in ichthyoplankton identification following the morphological feature collection. The development of a morphological database combining COI barcodes library of ichthyoplankton is needed for future investigation reference (Hubert et al., 2015;Hou et al., 2021a). In this study, we first preserved the ichthyoplankton samples in 5% neutral formalin and then transferred to 95% ethanol in less than 12 h. Highquality COI gene fragments obtained from eggs and larvae were 54.1 and 60.9%, respectively.
The 60 ichthyoplankton taxa identified in this study inferred the underestimation of the ichthyoplankton diversity in Sansha Bay because many fish larvae have a vertical distribution pattern on a diel basis (Ahlstrom, 1959;Rodríguez et al., 2006;Auth et al., 2007). Although we realized the importance of both day and night sampling for understanding the ichthyoplankton composition and assemblage, night sampling was not applied mainly for safety considerations. Moreover, adhesive eggs are not able to collect using plankton nets. In the future, sampling methods need to be modified to obtain more ichthyoplankton samples.
The abundance of ichthyoplankton species (proportion of individuals) in Sansha Bay has changed over time. For example, a high abundance was found in the black sea bream Acanthopagrus schlegelii (Bleeker, 1854) (22%) and the bighead hairtail  (Dai, 2006;Wang et al., 2010;Shen, 2011). The shift of major ichthyoplankton composition may be driven by intensive fishery pressure, habitats loss, or climate change (Bui et al., 2010;Fodrie et al., 2010;Shen, 2011;McCain et al., 2016;Curran et al., 2021). In contrast, the sampling frequency (monthly  or seasonally), spawning ground, and station setting may also influence the ichthyoplankton composition. The temporal and spatial variation of high ichthyoplankton density was detected from 1990 to 2019 (Figure 8) in Sansha Bay.

Spawning Activities in Sansha Bay
This study confirmed that Sansha Bay is an important area for fish reproduction, for pelagic, demersal, and benthic species, including some commercially important species. One of the most traditionally important fishery species in Sansha Bay was L. crocea, a well-known species from spawning and overwintering aggregations in Chinese coastal waters (Liu and Sadovy de Mitcheson, 2008). It was documented that L. crocea migrated into Sansha Bay in May and June to spawn through the narrow Dongchong Channel (Chu and Wu, 1985). Eggs and larvae of L. crocea were collected in May, June, August, October, and November (Dai, 2006;Shen, 2011;Xu, 2018), confirming the existence of two spawning seasons (spring and autumn) in the region (Liu and Sadovy de Mitcheson, 2008). In this study, eggs and larvae of L. crocea were collected in April−June because only spring spawning season was focused (Supplementary  Tables 1, 2), and mature females and mature males were also Explains %: the percentage of the total variation explained by the explanatory variable; contribution %: the percentage of the variable contribute to the explanatory power. caught in April−June by fixed nets in Sansha Bay following the confirmation by gonadal histology (Liu et al., 2020a). It is the first record that L. crocea also spawned in April. Further understanding the reproductive dynamic and migration pattern for L. crocea in the Sansha Bay will help us assess the stock status of this critically endangered species (Liu et al., 2020b). Two sciaenid species, namely, the bighead croaker Collichthys lucidus (Richardson, 1844) and N. albiflora, have been of commercial importance in Fujian waters (Chu and Wu, 1985). Juveniles and adults can be caught by fixed nets in Sansha Bay throughout the sampling period of April−September (Liu et al., 2020a). Ichthyoplankton samples (eggs as majority) of C. lucidus and N. albiflora were collected in May, June, and August and in April−September, respectively, with dominant in May and June, and in May−September, respectively (Supplementary Table 1 and Table 2). However, ichthyoplankton of C. lucidus and N. albiflora were not collected in previous studies in , 2007in (Dai, 2006Wang et al., 2010;Shen, 2011;Xu, 2018). Based on the external morphology identification of these previous studies, it was possible that the eggs and larvae of the two species were grouped as Sciaenidae species. In contrast, the reproductive dynamics of C. lucidus and N. albiflora have not been examined in Sansha Bay; therefore, it is unclear if their reproductive behavior has changed over time. For C. lucidus, the spawning and peak spawning seasons in Pearl River Estuary, approximately 800 km south of Sansha Bay, were December-July and May, respectively (Ou et al., 2012). For N. albiflora, the spawning and peak spawning seasons in Xiangshan Bay (Zhejiang Province), approximately 400 km north of Sansha Bay, were May−July and June-July, respectively (Lin et al., 2013). The variation on spawning seasons of the same species in different geological locations merits further assessments.
In Fujian waters, seabreams (Sparidae) such as A. schlegelii and the red seabream Pagrus major (Temminck and Schlegel, 1843) mainly spawned in February−April and October−December, respectively (Chu, 1985). The eggs and larvae of A. schlegelii were collected from April to August in Sansha Bay (Wang et al., 2010;Shen, 2011;Xu, 2018; this study), indicating the species has a longer spawning season, extending from spring to summer. For P. major, the eggs and larvae were collected in November 1990 (Dai, 2006) and in April and May 2019 (this study) in Sansha Bay, indicating the existence of spring and autumn spawning seasons of the species.
In the East China Sea, the spawning season of T. lepturus is March−August, producing pelagic eggs, and the spawning grounds are known offshore of Chinese waters (Chu, 1985;Nakamura and Parin, 1993;Xu and Chen, 2015). However, T. lepturus was dominant in ichthyoplankton samples (eggs) in May and August 1990 in Sansha Bay (Dai, 2006), inferring its possible spawning grounds in nearshore and semi-closed bays. It was not collected in 2007,2008,2010, and 2019 surveys covering its spawning season in Sansha Bay (Wang et al., 2010;Shen, 2011;Xu, 2018;this study). Nearly 100% of T. lepturus collected in Sansha Bay using fixed nets were young juveniles (3.2−46.1 cm total length, TL) in May−September 2019 (Liu et al., 2020a). Although the spawning ground function of T. lepturus in the Sansha Bay is still not clear, its nursery ground function is clear.
The Indo-Pacific king mackerel Scomberomorus guttatus (Bloch and Schneider, 1801) and the Japanese Spanish mackerel S. niphonius are pelagic, forming spawning migrations from deep waters to shallow waters in spring and spawning in May and June in Sansha Bay (Wu, 1985). Ichthyoplankton of S. guttatus and S. niphonius were collected in June 2010 (Xu, 2018), and in April and May 2019 (this study), confirming the spawning ground function of the Sansha Bay for the two species.
Ichthyoplankton of two alien species were first collected in Sansha Bay in this study; eggs of S. ocellatus in August    Table 1). The two species were introduced into China for mariculture purposes from United States in 1991 and from France in 2001, respectively, and juveniles and adults have dispersed in coastal waters of China (Wang and Ji, 1996;Liu et al., 1998;Wang et al., 2006;Lin et al., 2020;Zhang, 2020). It is likely that S. ocellatus has established reproductive wild stocks in Chinese waters (Lin et al., 2020). Both species are commonly cultured in floating cages in Sansha Bay, and their juveniles (S. ocellatus: 35.8 cm TL, n = 1; S. aurata: 3.7−31.2 cm TL, n = 9) were collected using the fixed nets in 2019 (Liu et al., 2020a). This study revealed that S. ocellatus and S. aurata can reproduce and may be able to establish spawning population in Sansha Bay.

Ichthyoplankton Assemblages
Five main temporal ichthyoplankton assemblages were apparently separately in Sansha Bay. The densities of ichthyoplankton along Dongchong Peninsula coastline into Dongwuyang (S12-S25) were higher than those in the northwest waters of Sandu Island (S01-S11) with the significant high salinity (this study) and high tidal currents speed (Lin et al., 2018), and the dominant species were from Engraulidae [S. tenuifilis], Soleidae (Solea ovata Richardson, 1846), Serranidae [Epinephelus akaara (Temminck and Schlegel, 1842)], Sciaenidae (N. albiflora and C. lucidus), and Sparidae (A. schlegelii and P. major), including pelagic, demersal, and benthic species. The results indicated that these species may prefer spawning in waters with high salinity and strong currents, or fish eggs and larvae were dispersed to there by currents.
The spatial distribution patterns of ichthyoplankton varied over the past three decades in Sansha Bay, and three zones can be considered as core areas (Dai, 2006;Wang et al., 2010;Shen, 2011;Xu, 2018; this study; Figure 8). First, it is around Sandu Island. The high ichthyoplankton density area was in the north in 1990, then, shifted to the east in the late 2000s, and was not important in 2019. Second, it is in Guanjingyang waters with an extension to the Dongchong mouth. The high ichthyoplankton density was found in 2008, 2010, and 2019. Third, it is in Dongwuyang waters. The high ichthyoplankton density was found in 1990, 2007, 2010, and 2019. Overall, the high ichthyoplankton density areas remain the same in Guanjingyang and Dongwuyang over the past three decades and reduce largely from 1990 in Sandu Island water.

Influence of Abiotic and Biotic Factors
Studies revealed that both abiotic and biotic factors shaped the fish assemblages (including fish reproduction stocks and ichthyoplankton) (Olivar, 1990;Marshall and Elliott, 1998;Pombo et al., 2005;Zhang et al., 2015;Hsieh et al., 2016;Rodriguez, 2019). Temperature and salinity influence fish distribution through the different preferences and tolerances of the species. Temperature is likely to be the dominant factor influencing the variability of migration, spawning, and recruitment (Gibson et al., 1993;Marshall and Elliott, 1998) and the abundance and distribution of phytoplankton, zooplankton, and ichthyoplankton (Siokou-Frangou et al., 1998;Rautio, 2001;Hsieh et al., 2005;Gogoi et al., 2021). In this study, seven abiotic and biotic factors, including T, S, DO, pH, Chla, Phy, and Zoo, significantly structured the ichthyoplankton assemblages in Sansha Bay, and temperature was also the dominant factor and then salinity and Chla.
The dominant species of phytoplankton and zooplankton in May and June (peak spawning season) in Sansha Bay were mainly from Bacillariophyta and copepods (Liu et al., 2020a), all of which may play important roles as feeds for fish larvae and young juveniles to further impact the distribution pattern of ichthyoplankton in Sansha Bay. Chla, representing the primary productivity, is associated with phytoplankton and zooplankton abundance and always has a positive correlation with phytoplankton abundance (Felip and Catalan, 2000). In this study, the red tide forming species Noctiluca scintillans (Macartney) (Kofoid and Swezy, 1921) bloomed in May 2019, which was without chloroplast and purely phagotrophic (Balch and Haxo, 1984;Zhang et al., 2020), resulting in apparently high abundance of phytoplankton and low concentration of Chla (Nakamura, 1998). Moreover, N. scintillas can feed on phytoplankton, zooplankton, and even fish eggs and larvae (Nakamura, 1998;Quevedo et al., 1999;Umani et al., 2004); therefore, the impact of biotic factor Chla related with N. scintillas might have a negative correlation with Phy, Zoo, and certain ichthyoplankton species (Figure 7).

CONCLUSION
In this study, 60 taxa ichthyoplankton were identified from horizontal and vertical tows by combining DNA barcoding and external morphology from April to September 2019 in Sansha Bay of Fujian Province, China. Among these, 52 taxa were identified at species level, 5 taxa at genus, and 3 taxa at family level. The 10 dominant species were from Engraulidae (S. tenuifilis and Stolephorus commersonnii), Serranidae (E. akaara), Leiognathidae (Photopectoralis bindus and Nuchequula nuchalis), Sciaenidae (N. albiflora and C. lucidus), and Sparidae (P. major and A. schlegelii). The temporal distribution of ichthyoplankton were divided into five assemblages, namely, April, May, June, July, and August−September with a high density occurred mainly in Guanjingyang and along the Dongchong Peninsula coastline into Dongwuyang. The temporal and spatial patterns were closely related to both abiotic and biotic factors. Temperature was the main factor influencing the assemblage structure and then Chla and salinity. Other factors, such as ocean currents, substrates, and other anthropogenic activities still need further investigation. Currently, little is explored to understand the entire life cycles of the dominant and commercially important fish species in Sansha Bay and its adjacent waters, the East China Sea. To achieve effective fishery stock management, understanding different life stages and their habitat uses, and the connectivity among different life stages are highly recommended.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: OK012006-OK012065).  The funders had no role in data collection and analysis, decision to publish, or preparation of the manuscript.