Assemblage Structure of the Ichthyoplankton and Its Relationship With Environmental Factors in Spring and Autumn off the Pearl River Estuary

Ichthyoplankton assemblages and their relationship with environmental variables are investigated in waters off the Pearl River Estuary in spring and autumn of 2019. Of 80 ichthyoplankton taxa identified using DNA barcode and morphological methods, 61 are identified to species. The most abundance families (Carangidae, Trichiuridae, Mullidae, and Scombridae) account for 61.34% of the horizontal total catch in spring, while Menidae and Carangidae are the most abundant families identified in autumn, accounting for 89.72% of the horizontal total catch. Cluster analysis identifies three species assemblages in spring, and four in autumn based on horizontal trawls. Relationships between assemblage structure and environmental variables (in situ and remote sensed) are determined by canonical correspondence analysis. Ichthyoplankton assemblage structure appears to be strongly influenced by sea level anomalies, salinity, water depth, temperature at 10 m depth, and distance from shore. We demonstrate the efficacy of using DNA barcode to identify ichthyoplankton, and suggest how these data can be used to protect fish spawning grounds in waters off the Pearl River Estuary.


INTRODUCTION
Information on the composition, and spatial and temporal distributions of ichthyoplankton assemblages is of value because it can be used to identify fish spawning sites and reproductive seasons, and appropriate areas and strategies for their protection (Potter et al., 1990;Raynie and Shaw, 1994;Oliveira and Ferreira, 2010;Makrakis, 2019). Ichthyoplankton-the early life stages of fish-undergo dynamic changes in morphology, physiological capability, behavior, ecological role, and habitat as they develop. Their natural mortality rate is high (Fuiman, 2003). Being a key life cycle stage, the survival rate of these fish eggs and larvae affects their recruitment and assemblage dynamics (Chambers and Trippel, 1997;Cao et al., 2007). The formation of ichthyoplankton assemblages is regulated by the physical and biological processes, and their distributions are considered to be more important to the survival of larvae than absolute larval abundance (Hewitt, 1981). Because ichthyoplankton assemblage distribution can be correlated with adult spawning phases and biotic/abiotic variables, those conditions suitable for adult reproduction, and mechanisms that affect planktonic early life stages, can be identified (Doyle et al., 1993;Whitfield, 1999;Hare et al., 2001;Mesa et al., 2016).
Ranked as the 13th or 14th largest river in the world and 2nd largest river in China, the Pearl River, has an annual runoff exceeding 33.6 × 10 9 m 3 (Ou et al., 2007). This river discharges into the South China Sea (SCS) across the Pearl River Delta floodplain, and in doing so transports abundant nutrients into the sea (Lu et al., 2007;Chen et al., 2008). The terrain of the Pearl River Estuary is complex, and the interlaced river network empties into the northern SCS through eight distributaries. This geographic complexity, coupled with variable salinity and high nutrient input, render the Pearl River Estuary and adjacent aquatic environments variable, changeable, and high productive (Zhang et al., 2011;Zhu et al., 2019). Waters of and off the Pearl River Estuary provide habitat for 542 species of fish (287 species in the estuary area, and 330 species in coastal and offshore areas) (Li, 2008). This species richness is more than 36% of the total number of fish species (nearly 1500) reported from the northern SCS (Zhan, 1998;Li, 2008;Sun and Chen, 2013). This high richness and the fish yield from these waters from depths greater than 30 m render this area one of the most important fishing grounds in the SCS (Wang, 2012;Cai et al., 2018). However, fish resources in this region have decreased dramatically in recent decades, largely because of anthropogenic disturbance and climate change (Zeng et al., 2019;Zhou et al., 2019). Meanwhile, the Indo-Pacific humpback dolphin, Sousa Chinese Osbeck, 1765, Indo-Pacific finless porpoise Neophocaena phocaenoides (Cuvier, 1829), Eden's Whale, Balaenoptera edeni Anderson, 1879, Bryde's whale, B. edeni Olsen, 1913, and humpback whale Megaptera novaeangliae Borowski, 1781, which mainly prey on fishes, have all been reported from waters off the Pearl River Estuary (Jefferson et al., 2012;Peilie, 2012;Lin et al., 2019). Sadly Bryde's and humpback whales have disappeared from these waters, and the population of Indo-Pacific humpback dolphin has continuously declined (Hung et al., 2006;Peilie, 2012;Lin et al., 2019). A shortage of food is considered to be one of the main reasons for these disappearances and population declines (Li et al., 2020).
Understanding the spatial distributions of fish is necessary to understand marine ecosystem dynamics. The Pearl River Estuary and its offshore waters are important spawning, feeding, and nursing ground in the northern SCS (Zhan, 1998;Li et al., 2000;Wang and Lin, 2006;Xiao et al., 2013). Accordingly, understanding how the distributions of ichthyoplankton are linked to environmental variables can provide valuable information to assist with spatial, temporal, and/or ecosystem-based approaches to fishery management Laman et al., 2017). However, despite the importance of these waters to fish stocks, due to the difficulties identifying early life stages of fish, only Xiao et al. (2013) has studied the ichthyoplankton fauna from them, wherein the abundance of 113 identified taxa was reported. Anthropogenic disturbance, climate change, and fisheries over-exploitation have all contributed to serious declines in fishery resources in the waters of and off the Pearl River Estuary (Duan et al., 2009;Wang et al., 2016;Zeng et al., 2019;Zhou et al., 2019). Because of this there is a pressing need to investigate the status of fish spawning grounds, and ichthyoplankton assemblages in this area, so that more effective management and/or protection of these resources can be achieved.
Precise identification of species in ichthyoplankton samples is important to understand species reproductive biology, preferred spawning time, and spawning grounds (Baumgartner et al., 2004;Cao et al., 2007;Claydon et al., 2014;Zhang et al., 2015). However, accurately identifying fish eggs and larvae to species is difficult, probably for which reason there have been few empirical studies on ichthyoplankton in the SCS (Li K. Z. et al., 2014;Huang et al., 2017;Hou et al., 2021). DNA barcode, which delimits species using molecular maker (e.g., cytochrome c oxidase subunit I, COI), have been widely applied to facilitate the identification of fish eggs and larvae (Valdez-Moreno et al., 2010;Frantine-Silva et al., 2015;Hubert et al., 2015;Leyva-Cruz et al., 2016;Burrows et al., 2018;Kerr et al., 2020;Chen et al., 2021;Hou et al., 2021). For instance, Leyva-Cruz et al. (2016) identified 42 taxa of fish eggs collected from the waters surrounding Banco Chinchorro in the Mexican Caribbean. Burrows et al. (2018) applied the DNA barcode to identify fish eggs in the Gulf of Mexico, and delimited 709 eggs to 62 species. Kerr et al. (2020) successfully identified 564 fish eggs to 89 taxa collected from northwestern Cuba and across the Florida Straits using DNA barcode. Therefore, it has been demonstrated that DNA barcode can effectively identify fish eggs and larvae, and illustrated the potentiality in the study of ichthyoplankton assemblage and ecology.
To identify fish spawning grounds in waters off the Pearl River Estuary we undertook ichthyoplankton surveys in the spring and autumn of 2019. Our aims were to: (i) identify fish eggs and larvae using DNA barcodes and morphology, and develop a morphological database for the ichthyoplankton in this area, (ii) describe the ichthyoplankton assemblages structure in spring and autumn, and (iii) reveal relationships between environmental variables and ichthyoplankton assemblage structure. By evaluating the influence of human activities on recruitment of fishery resources, the management and conservation of these species and their spawning grounds in waters off the Pearl River Estuary can be improved.

Study Area and Sample Collection
Surveys occurred at 12 stations in waters off the Pearl River Estuary (112.50-114.50 • E, 20.00-21.50 • N) in spring (April) and autumn (October) of 2019 (Figure 1). Ichthyoplankton samples were collected by zooplankton net (0.8 m diameter, 2.7 m long, 505 µm mesh, with cod-end container mesh of 400 µm) at stations over seabed depths exceeding 30 m. Four zooplankton net deployments were made at each station: two simultaneous horizontal tows at fixed depth (15 min at 1.5-2.2 knots), and two vertical hauls (∼sea bed to the surface, hauled at 1.5 m/s). To preserve useful morphological characteristics of fish eggs and larval specimens, and to ensure that we could extract quality COI sequences from samples, we (upon collection) immediately preserved ichthyoplankton samples in 4% neutral formalin for 2 h, then filtered this out and transferred fixed samples to ∼75% ethanol.
In the waters off the PER, several oceanographic events occurred, i.e., salinity front, coastal current, river plume, and water masses etc. (Dong et al., 2004;Su, 2004;Li et al., 2018), and these events were evaluated for possible effects on the ichthyoplankton assemblages. Temperature and salinity were considered as key variables in the oceanographic events. In the present study, the surveyed water was at the sea side off the salinity front (>32 PSU, the offshore boundary of the river plume; Ou et al., 2007). Thus, environmental variables included surface sea water temperature (SST, • C), sea water temperature at 10 m (T-10, • C), sea water temperature at 40 m (T-40, • C), sea surface salinity (SSS, PSU) were considered in the present study. In addition, sea level anomaly (SLA, mm), water depth (m), and offshore distance (km) were also adopted (He et al., 2014. SST and water depth were measured in situ. Satellite remote-sensed T-10, T-40, and SSS and SLA data were downloaded from HYCOM model data 1 , for which temporal resolution was 1 days, and spatial resolution was 1 km. SLA data were derived from sea surface elevation data of HYCOM, then calculated using the Google Earth Engine, a cloud platform by mean processing algorithm. Minimum offshore distance was calculated in R version 4.0.5 (R Core Team, 2021) using the Sp package (Pebesma and Bivand, 2005).

Laboratory Analyses
Ichthyoplankton were identified in the laboratory by combined the molecular and morphological analyses. Briefly, ichthyoplankton specimens were first identified by cytochrome c oxidase subunit I (COI) to the lowest taxonomic level using a combined local DNA barcode library (Hou et al., 2018) and the Barcode of Life Data (BOLD) system, and the blast search criteria similar to those used in Hubert et al. (2015). Because formalin-fixation can reduce the success rate of DNA extraction and sequence amplification after 10 days (Hou et al., 2021), especially for fish eggs, we quickly picked out the ichthyoplankton specimen beneath a stereomicroscope and randomly selected 15 eggs/larvae in each horizontal and vertical trawl to be DNA extracted and amplified for the first time. If possible we performed this 2 or 3 times for each sample if a sample contained >60/90 specimens. In the 3 times, the success rate of DNA extraction and sequence amplification reduced sharply to below 20%, after which we stopped. For the purpose of developing a morphological database for the ichthyoplankton in the study area, each specimen was photographed by a Zeiss microscope (Axioplan 2 imaging E) before the DNA extraction (Figure 2). For the specimens that failed in the DNA extraction were identified by morphology following Okiyama (1989) and Shao et al. (2001).

Data Analysis
Species density (abundance per tow) was standardized to catch per unit effort. Taxa were assigned to different ecological guilds using habitat data for adult fish species in Fishbase 2 . An index of relative importance (IRI) was calculated for each species or taxon using the following function (Zhu et al., 2002;Ren et al., 2016): IRI = N% × F% × 10,000 Where N% and F% are the relative abundance of the total catch and frequency of occurrence. Accordingly, each species or taxon was then defined as dominant (IRI ≥ 100), common (10 ≤ IRI < 100), or rare (IRI < 10).
Species accumulation curves were used to assess if our sampling effort was sufficient to describe ichthyoplankton species richness. The classic method "random" was used as an accumulator function to determine the means and standard deviations of species accumulation curves from random subsampling of the data without replacement (Gotelli and Colwell, 2001). The species accumulation curves were fitted in the "vegan" package implemented in R version 4.0.5 (Oksanen et al., 2020;R Core Team, 2021).
Species/taxa which accounted for more than 1% of the total catch in averaged horizontal tows were used in seasonal assemblage structure analysis. Data were first square root transformed and standardized to enhance the weighting of rare species (Clarke, 1993), and then clustered using Bray-Curtis similarity (Primer Software, Version 5.2). A non-parametric ANOSIM analysis was performed to examine differences within clusters in seasons (Clarke and Warwick, 2001). SIMPER analysis was used to calculate the contribution of each taxon to the Bray-Curtis similarity within clusters, and dissimilarity among clusters (Clarke and Warwick, 2001). Canonical correspondence analysis (CCA) was used to elucidate relationships between environmental variables and ichthyoplankton assemblages (CCA, CANOCO Software, Version 4.5). A category of sample coded "no fish" was created to prevent CANOCO from eliminating samples containing no ichthyoplankton (Grothues and Cowen, 1999). Down-weighting of rare species was performed. A forward selection was applied to test the statistical significance of environmental variables by a Monte Carlo permutation test (999 permutations) (Ye et al., 2015). Variation inflation factor (VIF) was used to test for collinearity between independent environment variables, with an upper cutoff of 10 (indicating excessive multicollinearity) (Graham, 2003;Kutner et al., 2004). Ordination diagrams of the CCA results were plotted to demonstrate associations between environmental variables and species. The length of a vector of a given variable on the CCA plot indicates the importance of that variable, with any variable having a correlation coefficient ≥ |0.4| is considered to be biologically important (Rakocinski et al., 1996;Akin et al., 2005). Species plotted closer to any vector with a longer distance from the origin have a stronger relationship with that vector (Ren et al., 2016). Species located near the origin either do not show a strong relationship to any of the variables, or they are found at average values of environmental variables (Elliott, 1998;Akin et al., 2005).

Ichthyoplankton Species Identification by DNA Barcode and Morphology
Of 4769 ichthyoplankton specimens collected, 2867 were from spring and 1902 from autumn. In all, 931 fish eggs and 229 larvae were selected for photography and DNA extraction, from which high-quality COI sequences were obtained from 308 fish eggs (33.1% of samples) and 125 larvae (54.6% of samples). A total of 75 taxa were identified using BOLD and local DNA barcode libraries (Supplementary Table 1). Among the 75 taxa identified by COI sequences, 61 were identified to species (referred to 50 genera, 33 families, and 10 orders), 5 were identified to genus, and 9 could not be identified. Of these, 38 taxa identified from eggs were attributed to species (comprising 31 genera, 20 families, and 8 orders), and 35 larvae were identified to species (comprising 29 genera, 21 families, and 7 orders) (Supplementary Table 1). In all, 80 taxa were identified by morphological and molecular analyses (75 taxa by COI sequences, 4 by morphology, and 1 taxon that are unidentified) (Supplementary Tables 1, 2). Our species accumulation curves were non-asymptotic, indicating substantially more taxa would have been encountered with increased sampling effort in spring and autumn (Figure 3).

Spatial and Temporal Patterns of Ichthyoplankton in Spring and Autumn
In spring, in horizontal net tows, 11 taxa were categorized as dominant, contributing 82.46% to the total abundance in these samples. In vertically hauled net tows, nine taxa were dominant, contributing 77.18% to vertical samples. In autumn, in horizontal net tows, five taxa were dominant, contributing 92.92% to the total abundance of species in these samples. In vertically hauled net samples, three species were dominant, contributing 79.52% to these samples (Table 1 and  Supplementary Table 2). Among these taxa, Terapon jarbua (Forsskål, 1775) was dominant in both horizontal and vertical trawls during spring and autumn.
The mean total abundance in spring was 91.7 ind./net in horizontal tows and 27.75 ind./net in vertical hauls. These values were higher than those in autumn: 68.88 ind./net (horizontal tows) and 10.38 ind./net (vertical hauls). However, One-way ANOVA indicated that there are no significant difference for the abundance in horizontal and vertical tows between spring and autumn (F spring = 0.19, F autumn = 2.08, P > 0.05). In the present study, we took horizontal samples to examine differences within clusters in seasons. In spring, ichthyoplankton composition was not very similar, and there was significant    Table 3). In spring, there were three main clusters of taxa comprising 18 categories of ichthyoplankton, each represented more than 1% of the total abundance (I to III, Figures 4A,B). In autumn, there were four clusters indicated by six categories of ichthyoplankton (Figures 4C,D).

Relationships of Ichthyoplankton Assemblages to Environmental Variables
Environmental variables off the PER are shown in Supplementary Table 4. Because a VIF test indicated that SST and T-40 showed multicollinearity with other environmental variables, we excluded them from subsequent analysis, leaving us with five variables (SLA, SSS, T-10, Depth, and Offshore distance) for CCA analysis. Monte Carlo permutation tests indicated that four (excepting SLA) of these variables contributed significantly to explaining ichthyoplankton assemblage structure (P < 0.05). According to the rule of correlation coefficient ≥ |0.4|, SLA was evaluated to be biologically important and used in CCA analysis (Supplementary Table 5). Eigenvalues were 0.413 (CCA1), 0.295 (CCA2), 0.132 (CCA3), and 0.095 (CCA4). In brief, the first and second CCA axes accounted for 41.6 and 29.8% of the variance in the species-environment relationship, respectively, with the correlation efficiencies of the first two axes being 0.878 and 0.921, respectively ( Table 3). The sum of all canonical eigenvalues (0.992) only equaled 33.09% of the unconstrained eigenvalues (2.998), showing the restrictive effect of building environmental relationships into the CCA model. Results obtained from the first two axes were plotted to explain species-environment relationships. The first axis was highly correlated with T-10 and distinguished autumn sampling stations. The second axis was highly correlated with SSS, SLA, and water depth, and distinguished spring sampling stations ( Figure 5A). The ordination diagram of the first two axes with the scores for environmental variables and species showed that most species scattered near the origin, which indicated average values in relation to environmental variables for these species. However, some species were outliers, e.g., Pristigenys niphonia (Cuvier, 1829), Istiophorus platypterus (Shaw, 1792), and Elagatis bipinnulata (Quoy and Gaimard, 1825), which was separated from the majority of species along the first axis, showing mainly a negative correlation with SSS; Euthynnus affinis (Cantor, 1849), Ophisurus serpens (Linnaeus, 1758), Tentoriceps cristatus (Klunzinger, 1884), and Unidentified Gobiidae, which were separated from most species along the second axis, showing mainly a positive correlation with SLA and a negative correlation with water depth and offshore distance ( Figure 5B).

Performance of Molecular Techniques for Identifying Ichthyoplankton
We identified 75 taxa using DNA barcodes, combining the BOLD system with a local DNA barcode library. Of these taxa, 61 were identified to species, including 38 taxa from eggs, and 35 taxa from larvae; we could not attribute 11 egg or 6 larval COI sequences (3.93% of all sequences) to any taxon, indicating incompleteness of regional barcode libraries (Supplementary Table 1). Our study demonstrates the importance of DNA barcoding in studies involving ichthyoplankton identification, especially their egg stages. However, for sequencing large numbers of ichthyoplankton samples this technology is expensive and time consuming, and its application would be impracticable for large spawning-ground surveys. To obtain species-level information from ichthyoplankton using molecular approaches, it proved indispensable to first differentiate morphotypes based on their morphology. Thus, developing a morphological database together with a corresponding DNA barcode library will facilitate identification of fish early life stages (Hubert et al., 2015;Hou et al., 2021). In the present study, we develop a preliminary morphological database of 75 ichthyoplankton taxa for these waters, layed the foundation for the morphological identification in the waters off the PER (Figure 2). Compared to the previous study, we improved the preserved method, transferred from the 4% neutral formalin to ∼75% ethanol, which can ensure that we can obtain a certain proportion of COI sequences of both fish larvae and fish eggs, although the morphology quality of larvae is not good as those preserved in formalin solution (Hou et al., 2021). In addition, the present study also represents the first attempt by using the DNA Barcode to provide more detailed information regarding the species composition, occurrence of spawning site and spawning seasons for ichthyoplankton off the PER. However, the present surveys effort were not robust estimates of species richness due to the lack of asymptote in the species accumulation curves (Figure 3), indicating that adding more survey stations would added many more ichthyoplankton species.

Spawning Periods and Ichthyoplankton Assemblage Structure off the Pearl River Estuary
Because we surveyed only in spring (April) and autumn (October), we are limited to discussing these two sampling events. The 33 ichthyoplankton taxa that occur in both spring and autumn samples account for 41.25% of all species (80 taxa) that we report from this area. This suggests that many fish species spawn in both spring and autumn. Most species also have spawning peaks in spring, which corresponds to subtropical monsoon characteristics in the northern SCS and adjacent waters (Tzeng et al., 2002;Hsieh et al., 2005;Hsieh H. Y. et al., 2010;Xiao et al., 2013;Li K. Z. et al., 2014). Earlier analyses of fish spawning times were based largely on monitoring gonad development of fishes in the northern SCS (i.e., begin at gonadal maturity stage IV); such an approach is also expensive and requires considerable manpower and logistical resources (Chen et al., 2002. While spawning periods could be inferred from earlier monitoring work, the spawning locations of fishes remained unknown. Thus, that information we present from this region provides clear evidence of the spatial and temporal  distribution of ichthyoplankton in spring and autumn, and likely presence of spawning grounds nearby. To protect spawning fish stocks and juveniles, a 2-month (June and July) moratorium on fishing in the SCS was implemented from 1999 to 2008; and extended to 77 days from 16 May to 1 August from 2009 to 2016; since 2017 it was even extended to 3.5 months (1 May to 16 August). Our study indicates that the moratorium in effect from 2017 does not protect spawning fish in this region in spring. For 1 month prior to closure, fishers off the Pearl River Estuary could be catching spawning adults, which intuitively would affect recruitment and exacerbate resource declines. Within Pearl River, a fishing moratorium extends from 1 March to 30 June-a closure which coincides with the spawning periods of freshwater fish, effectively increasing their recruitment (e.g., Megalobrama hoffmanni Herre and Myers, 1931) (Li Y. F. et al., 2014). Our data for offshore waters support advancing the fisheries moratorium in waters off the Pearl River Estuary by at least 1 month to April to better protect fishery resources.
Because we used a zooplankton net to collect drifting and suspended eggs, many eggs (e.g., demersal eggs of sea catfish, mouth-hatched eggs such as those of cardinalfish, and adhesive eggs) were not sampled. Accordingly the composition of our ichthyoplankton assemblages is somewhat biased, and not representative of the entire fish community known from these waters. However, we do report several unexpected and most important taxa, such as commercially important Scombridae species [Thunnus tonggol (Bleeker, 1851), E. affinis (Cantor, 1849), Auxis thazard (Lacepède, 1800), Auxis rochei (Risso, 1810)], and Nemipteridae species [Nemipterus virgatus (Houttuyn, 1782) andNemipterus bathybius Snyder, 1911], the most caught fish taxon Decapterus spp. in the northern SCS, and the new species Strophidon tetraporus sp. nov. (Huang et al., 2020). The spawning grounds of these species are not well known, mainly because of difficulties experienced with accurate identification of their eggs and larvae, and lack of any illustrated handbooks. Additionally, benthopelagic fish (e.g., Diaphus watasei Jordan and Starks, 1904, and Cyclothone sp.) were identified from COI sequences, indicating that waters off the Pearl River Estuary are affected by continental shelf surface water mixing (Zhang et al., 2011).

The Environmental Factors That Affected the Ichthyoplankton Assemblages
Understanding relationships between ichthyoplankton assemblages and physical and biological processes is becoming increasingly important for ecosystem-based fisheries management. Physical process such as the estuaries (Zhang et al., 2015), ocean currents (Doyle et al., 2002;Sassa et al., 2007;Mullaney et al., 2011;Thompson et al., 2014), water mass (Grothues and Cowen, 1999;Hare et al., 2001;Espinosa-Fuentes and Flores-Coto, 2004;Muhling et al., 2010), monsoons (Hsieh et al., 2005;Hsieh H. Y. et al., 2010) can influence the distributions and survival of fish eggs and larvae, and define the structure and diversity of their assemblages. In the present study, the suite of environmental variables that we examined were significantly associated with ichthyoplankton assemblage structure. The first CCA axis (significantly influenced by SSS and T-10) represents a spatial gradient from the distance off the Pearl River Estuary. Increased salinity and temperature at 10 m depth were associated with changes in ichthyoplankton assemblage structure. Salinity and temperature affect fish reproductive activity and the spatial and temporal distribution of ichthyoplankton assemblages (Whitfield, 1999;Fuiman, 2003).
Salinity is inversely related to egg diameter through its influence on the width of the perivitelline space, and it also affects the hatching process and development rate (Fuiman, 2003). Temperature, a potent environmental regulator of fish physiology, can affect production and activity of hatching enzymes, and is strongly, negatively correlated with the duration of incubation (Fuiman, 2003). Thus, temperature and salinity drive reproductively mature fishes to appropriate habitat for spawning. In the nearshore and offshore waters of Pearl River Estuary, the salinity gradually increases from upstream to downstream, ranging (multi-year) 0.03-33.08 PSU through dry and wet seasons (Jia et al., 2011), meanwhile the sanility front discriminate the nearshore waters (with a salinity of <26, the nearshore boundary of the plume, Wong et al., 2003) and offshore waters (with a salinity of >32, the offshore boundary of the plume, Ou et al., 2007). The horizontal structure and shapes of Pearl River buoyant plume varied among monthly and annually, which are driven by the river discharge in different months and years (Ou et al., 2007;Li et al., 2018). This variability provides varied habitat for many (542) fish species (Wang and Lin, 2006;Zhu et al., 2019), which may also accordingly affect the spatial distribution and ichthyoplankton assemblages for this area. In the present study, the surveyed water was at the sea side off the salinity front (>32 PSU), the dominant characteristics taxa in clusters were associate with high salinities and temperature, i.e., Selar crumenophthalmus (Bloch, 1793), D. macrosoma Bleeker, 1851, A. thazard and A. rochei in spring, and M. maculata and Coryphaena hippurus Linnaeus, 1758 in autumn (Table 1 and  Supplementary Table 3).
Biological processes, i.e., spawning behavior of adults, duration of the incubation period, metamorphosis, interspecific competition, and predation may also affect ichthyoplankton assemblages (Harden Jones, 1968;Mason and Brandt, 1996;Miller, 2002;Olivar et al., 2010). Fish have evolved a wide range of life-history strategies in response to the environments, and the offspring must be produced at a specific time and locations, delivered to appropriate nursery grounds to survive and grow, and then mature to join the reproductive population at the appropriate place (Harden Jones, 1968). The interspecific competition for zooplankton prey may exist in the larval stages and lead to reduce growth and survival of competitors' larvae, which may disrupt the assemblages. Meanwhile, the predators can reduce the relative larval populations. In the nearshore and offshore waters of PRE, the abundant nutrients transported by river discharges produce high productivity of zooplankton (Li et al., 2006), which attract large number of fishes to spawn and nurse in the coastal waters, i.e., the D. maruadsi (South China Sea Fisheries Research Institute, 1966). As the D. maruadsi was the dominant prey for genus Trichiurus fishes (Yan et al., 2011), the prey availability induced the formation of predators' spawning ground, and then affected the ichthyoplankton assemblages.
Furthermore, variation in species distribution explained by the first four CCA axes was 31.2%, indicating that other factors influence the distribution and composition of ichthyoplankton assemblages off the Pearl River Estuary. The effect of other variables (e.g., chlorophyll-a, and phytoplankton and zooplankton concentration) on ichthyoplankton assemblage structure should be examined in subsequent investigations.
Our analysis of ichthyoplankton assemblages indicates that waters off the Pearl River Estuary are important spawning grounds for many commercially important species. We believe that consideration should be given to protecting this area.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The animal study was reviewed and approved by Guangdong Ocean University.

AUTHOR CONTRIBUTIONS
GH and HZ analyzed the data and completed the first draft. GH, YC, and JW performed the field survey. GH, JW, LL, YC, and JL conducted the laboratory experiments. CP and HZ provided the guidance on data analysis and structure. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank the editor and reviewers for their constructive comments on our manuscript. We also thank numerous members on the survey ship for their help.