Variation in Larval Fish Assemblage Dynamics Associated with Harmful Algal Blooms in a Temperate Estuary, South Africa

Eutrophication and concomitant harmful algal blooms (HABs) are on the rise globally and pose a threat to larval stages of fishes that rely on estuarine nursery habitats for growth and survival. The anthropogenically altered low inflow estuary (LIE), Sundays Estuary, South Africa, supports persistent and predictable HABs. This study investigated the effects of HABs on the environmental conditions and larval fish assemblage of this warm temperate nursery area. Sampling took place during the austral spring of 2016 and 2018 at three sites in the mesohaline zone where both larval fish abundance and HABs are known to reach maxima. Physico-chemical and phytoplankton data were collected during the day and night, while larval fishes were sampled after nightfall. Physico-chemical parameters and larval fish assemblages were compared by water column (surface and bottom) and sites within the mesohaline zone, years, and HAB intensity (hypereutrophic ≥ 2781 H. akashiwo cells•mL−1; standard ≥ 205 H. akashiwo cells•mL−1; none < 205 H. akashiwo cells•mL−1). A longer period of consecutive hypereutrophic bloom conditions was recorded during 2018 compared to 2016. Dissolved oxygen concentration was notably higher during hypereutrophic blooms (χ2 = 23.759, df = 2, P < 0.001) and reached a maximum of approximately 21 mg•L−1 during the day and 13 mg•L−1 at night. Density and similarity of estuarine resident larval assemblages were negatively correlated to supersaturated dissolved oxygen concentrations. Greater mean densities of estuarine resident larvae were recorded during hypereutrophic blooms compared to standard blooms and bloom absence and diversity was lower during 2018 when hypereutrophic bloom conditions were more persistent. These changes may have major implications for successful early development of fishes that rely on the Sundays Estuary and similar LIEs as a nursery.


Introduction
Habitat suitability of nursery areas is a critical requirement for larval fish survival and eventual recruitment and replenishment of adult breeding populations (Bailey and Houde 1989;Claramunt and Wahl 2000).Estuaries in South Africa provide ideal nursery habitats that support numerous fish species of conservation and economic importance due to higher productivity and more refuge opportunities than neighbouring marine waters (Whitfield 1994a;Strydom 2015).However, estuaries are often closely associated human settlements and are subject to deleterious upstream activities (Paerl et al. 2006), resulting in significant habitat degradation in estuaries worldwide (Van Niekerk et al. 2013;Strydom and Kisten 2020).By improving knowledge of these impacts, informed practices like improved water management, environmental impact assessments, restoration ecology, and conservation management can be implemented.This is particularly pertinent in a time where fish populations are facing collapse due to overfishing and other anthropogenic pressures.
Larval fishes require a suite of conditions to maximise their recruitment success in nursery areas.For example, temperature regulates metabolic processes and, in turn, growth rate (Houde 1989), while turbidity affects behaviour such as prey capture and predator avoidance (Hunter 1981;Lehtiniemi et al. 2005).Therefore, environmental disturbances that influence water quality and biotic conditions may in turn affect the growth and survival of fishes, although many of the effects are not wellunderstood.The prevalence of harmful algal blooms (HABs) across the globe is closely linked to eutrophication (Glibert 2017(Glibert , 2020)).HABs alter the water quality by means of oxygen depletion and saturation (Lemley et al. 2018b), exposing grazers to toxins or other harmful by-products (Almeda et al. 2011;Genitsaris et al. 2019) and reducing light availability to benthic habitats.There has been research on the effects of hypoxia and oxygen supersaturation on larval fishes, although these have mostly been limited to laboratory-based experiments (Breitburg 1994;Kolesar et al. 2010;Dong et al. 2011).The effects of toxins produced by phytoplankton are also rarely investigated in larval fishes, with a greater focus on adult fishes and invertebrates (Yu et al. 2010;Basti et al. 2016).
Low-inflow estuaries are particularly susceptible to anthropogenic influences (Warwick et al. 2018).On the south-east coast of South Africa, the low-inflow Sundays Estuary has received particular attention in terms of eutrophication.The estuary is situated in an arid region experiencing seasonal rainfall and is influenced by agricultural return flows and low, regulated river inflow resulting in predictable, stable conditions.These conditions have led to increasing HAB occurrence since the 1980s (Emmerson 1989;Hilmer and Bate 1990).In recent years, the recurrent and highly predictable nature of seasonal HABs have been documented (Lemley et al. 2017b).Besides eutrophication and HABs, the estuary has been the focus of several studies on zooplankton (Wooldridge and Bailey 1982;Jerling and Wooldridge 1995;Sutherland et al. 2013) and fishes, including larvae (Harrison and Whitfield 1990;Sutherland et al. 2012;Strydom et al. 2014) and juveniles (Beckley 1984).The predictability and well-studied nature of the estuary provided an ideal opportunity to investigate the effects of HABs on habitat suitability for larval fishes in a low-inflow estuary.
Heterosigma akashiwo, a flagellate belonging to the class Raphidophyceae, forms large blooms during the warmer months (spring/summer) in the Sundays Estuary (Lemley et al. 2018b).The species is mixotrophic (Zhang et al. 2006;Jeong 2011), tolerant of a broad range of salinities and temperatures (Martínez et al. 2010), and can migrate vertically (Bearon et al. 2004;Lemley et al. 2018a).It can also produce cysts which allow H. akashiwo to survive unfavourable conditions (Kim et al. 2015).This bloom-forming species is well-known in fish and shellfish aquaculture internationally for its adverse effects on cultured organisms (Wang et al. 2006;Yu et al. 2010).In association with H. akashiwo blooms in the Sundays Estuary, recent research reported poorer biochemical body condition of larval Gilchristella aestuaria (Smit et al. 2023), avoidance of blooms by adult Mugil cephalus (Bornman et al. 2021), gill alterations (Bornman et al. 2022a), and decreased abundance (Bornman et al. 2022b) in juvenile mugilids.
The aim of this study was to investigate the effects of HABs on the environmental conditions and larval fish assemblage of the Sundays Estuary.The specific objectives were to establish linkages between the environmental conditions of the water column (phytoplankton and physico-chemical) and the responses of larval fishes (density, diversity, and developmental composition and assemblage structure) under bloom conditions.It was predicted that larval density and diversity would be lower during extreme densities of the bloom species, H. akashiwo, associated with persistent bloom formation.

Study Area
The permanently open Sundays Estuary (33° 43′ S, 25° 51′ E) (Fig. 1) is situated along the warm-temperate, southeast coast of South Africa, and is characterised by a shallow (~ 2 m), channel-like system (Hilmer and Bate 1990), a stable, regulated hydrodynamic environment, and nutrientrich agricultural return flows (Emmerson 1989), resulting in a permanently eutrophic state (Lemley et al. 2017b).An inter-basin water transfer scheme has resulted in consistent low freshwater inflow throughout the year and a largely suppressed frequency of major natural pulse events (Lemley et al. 2017b).This consistent low input of freshwater, combined with semi-diurnal tides, maintains a full salinity gradient and permanent connection to the sea (Wooldridge and Bailey 1982).Vertical stratification and mesohaline conditions are also present year-round due to these modifications, while flushing events are intermittent due to flow regulation (Lemley et al. 2018a).The eutrophication of the estuary has led to seasonal and recurrent HABs which tend to be monospecific, often exceeding 80 µg•L −1 (Hilmer 1990;Hilmer and Bate 1990;Lemley et al. 2017bLemley et al. , 2018a)).

Field Methods
The mesohaline zone (5-18) in South African estuaries supports the plankton productivity maxima (Strydom et al. 2003(Strydom et al. , 2023)).The Sundays Estuary mesohaline zone is also the most conducive to HABs where H. akashiwo blooms prefer the lower salinity surface waters within the mesohaline zone (Lemley et al. 2017b(Lemley et al. , 2018a(Lemley et al. , 2021)).Therefore, variable bloom conditions are found within the mesohaline zone which coincide with the secondary production of the estuary (Strydom et al. 2003).Samples were collected at three sites within the mesohaline zone, allowing for a possible broader range of HAB biomass to be captured.To account for the tidal shifts during the study period, the study sites were determined by salinity (lower mesohaline (LM) ~ 5, middle mesohaline (MM) ~ 10, upper mesohaline (UM) ~ 18), rather than by predetermined fixed locations.Abiotic and biotic variables were assessed isochronously twice a week during austral spring in 2016 and 2018.Phytoplankton samples were collected during the day, when peak bloom activity occurred (Lemley et al. 2018a) alongside dissolved oxygen concentrations which were highlighted as an important variable during a pilot study (Smit et al. 2021).Nighttime sampling commenced 30 min after nightfall, during which physico-chemical data and phytoplankton and larval fish samples were collected at each site.Data collected during a pilot study from the MM site using the same sampling protocol in spring 2016 were supplemented to the 2018 dataset (Supplementary Material Table S1).

Physico-Chemical Parameters
Physico-chemical parameters were measured using a YSI Pro DSS 6600 multiparameter probe at 0.5 m intervals throughout the water column, i.e. surface to bottom-waters.Temperature (°C), salinity, and turbidity (NTU) were taken during the night while dissolved oxygen (mg•L −1 ) was recorded at day and night to capture exposure of the larval fishes to biologically stressful oxygen concentrations.These measurements were averaged at each site to match the section of the water column sampled for larval fishes, with the surface representing measurements taken at 0-0.8 m, while the bottom tow was represented by the depth of between 0.9 and 1.5 m depending on channel depth.

Phytoplankton Biomass
Phytoplankton biomass (µg chlorophyll-a•L −1 ) was determined using replicate water samples (n = 2) taken with weighted pop-bottles at the sub-surface (0 m), 0.5 m, 1.0 m, and bottom of the water column during the day.The water samples were gravity filtered in the field through Munktell MGC (1.2 µm pore-size) micro-glass fibre filters, which were subsequently frozen for lab analyses.Chlorophyll-a was extracted in 10 mL of 95% ethanol for 24 h in a cold (1-2 °C), dark room.After the 24-h extraction period, the extract was refiltered through Munktell MGC filters and absorbance of the filtrate was determined before and after acidification with 1 N HCl using a spectrophotometer (Nusch 1980).As with the physico-chemical variables, data were averaged according to depth (surface 0-0.8 m; bottom 0.9-1.5 m) at each site to match the water sampled by plankton nets for larval fishes.

Phytoplankton Enumeration and Community Composition
Samples were collected in the field using weighted popbottles at the sub-surface (0 m), 0.5 m, 1.0 m, and bottom of the water column and preserved with 1 mL of a 25% glutaraldehyde solution for laboratory analyses.Rose Bengal stain was added to a known volume of the preserved samples and were settled overnight in 26.5-mm-diameter Utermöhl settling chambers using the Coulon and Alexander (1972) method.Phytoplankton cells were counted and identified under an inverted Leica DMIL phase contrast microscope at 630 × magnification.A minimum of 200 cells or 200 frames were counted and used to determine the density (cells•mL −1 ) of species as per the equation described by Snow et al. (2000).The species were grouped into either diatoms or 'flagellates' to estimate the diatom:flagellate abundance ratio, which is an important consideration in terms of prey quality available to larval fishes (Black et al. 2016).

Larval Fish Enumeration and Community Composition
Two modified WP2 plankton nets (570 mm mouth diameter; 200 µm mesh aperture size) fitted with flowmeters (Kahlsico 005 WA 130 in 2016and General Oceanics 2030R in 2018) were used for plankton tows to collect larval fishes (Strydom 2015).The near-surface and near-bottom waters were sampled simultaneously on either side of the boat.The nets were lowered into the water simultaneously and towed for 3 min at 1-2 knots.The bottom net was held down using a graduated pole and typically ranged from 1 to 1.5 m in depth.Where possible, an oblique course across the channel of the estuary was followed to allow samples to be collected in the midchannel and near the margins (Strydom et al. 2002).After each tow, flowmeter readings were recorded to determine the volume of water filtered by the nets.Samples were fixed on site in 10% buffered formalin.
In the laboratory, larval fishes were identified and enumerated using a Leica M50 stereomicroscope.All larvae in samples were positively identified to the nearest taxonomic level using methods described by Neira et al. (1998).Larvae were classified into appropriate guilds according to the species dependence on estuaries (Potter et al. 2015).
The densities of larval fishes were calculated as follows and expressed as the number of individuals per 100 m 3 : where N = number individuals per sample, r = flowmeter revolutions, and c = predetermined calibration value (m 3 ) (Wooldridge and Erasmus 1980).

Phytoplankton
Phytoplankton were classified into three types of bloom presence ("none", "standard", and "hypereutrophic") based on Lemley et al. (2017a) which was subsequently used to compare larval fish assemblage and physico-chemical conditions.Blooms were defined according to their phytoplankton biomass (chl-a), while the corresponding density of the focus species, H. akashiwo, was used to indicate the HAB species presence.A regression analysis was performed between H. akashiwo cell density and log 10 -transformed chl-a concentration.The equation of the resultant trendline (3-order polynomial) was used to calculate the H. akashiwo density at "standard" (≥ 20 chl-a µg•L −1 ) and "hypereutrophic" (≥ 80 µg•L −1 ) blooms (Fig. 2) (Lemley et al. 2017b).Standard blooms were ≥ 205 H. akashiwo cells•mL −1 , and hypereutrophic blooms were ≥ 2781 H. akashiwo cells•mL −1 .Spearman rank correlation analyses Larval fish density no. per 100 m 3 = (N∕(r∕c)) × 100 were performed between the physico-chemical variables and H. akashiwo density to determine relationships.

Larval Fish Density
Larval fish density was grouped into estuarine resident and marine-spawned species according to Potter et al. (2015).Density of estuarine residents and marine-spawned species were compared by year, site (LM, MM, UM), water column (near-surface and near-bottom), and bloom presence using Mann-Whitney and Kruskal-Wallis tests.The relationship between larval fish density and environmental conditions was explored using generalised additive modelling (GAM).Initial models were built to include categorical factors, depth (surface and bottom) nested in sites (LM, MM, UM), and thin-plate regression splines of covariates (temperature, salinity, turbidity, dissolved oxygen concentration, H. akashiwo density, and diatom-flagellate ratio).Chl-a concentration was excluded since it was co-linear with H. akashiwo (P < 0.001; ρ = 0.654).Estuarine resident larval fish data followed a beta distribution which was not suitable for the GAM formula in mgcv in R; therefore, the data was log(x + 1)-transformed, which resulted in normal distribution, and subsequently, the GAM on estuarine resident density used the Gaussian family with log-link function.Marine-spawned larvae was zero-inflated, giving a negative binomial distribution.Data transformation on the marine-spawned larval data did not improve the distribution, and untransformed data were used in the model with the negative binomial family and log-link function.The dredge function from the MuMIN package in R was used to run possible models and produce the best models based on AIC (Akaike Information Criterion), whereby models with the lowest AIC were reported.This was performed in R statistical software (R Core Team 2019) using packages, MuMIn (Barton 2019) and mgcv (Wood 2011).Larval Fish Community Structure Diversity of estuarine resident and marine-spawned larval fishes was calculated using the Shannon-Wiener diversity index.Prior to analyses, larval fish density data were log(x + 1)-transformed.A Bray-Curtis similarity matrix by group-average was created for estuarine resident and marinespawned species.The matrix was used to perform and plot non-metric multi-dimensional scaling (NMDS).The correlation between environmental variables and group-differences were calculated using the function, envfit, in the R package, vegan, and plotted as overlays onto the NMDS plots, where only variables at a significance of P < 0.05 were included.Analysis of Similarity (ANOSIM) was performed to determine the differences between or among groups (year, site, water column, bloom presence).The contribution of species to these differences was calculated at a significance level of P < 0.05 using similarity of percentages (SIMPER).In terms of the analysis of the marine-spawned community, empty samples (samples containing no marine-spawned individuals) prevented the Bray-Curtis dissimilarity index from being calculated.Since empty samples were considered valuable and meaningful to this study, a zero-adjusted Bray-Curtis method was used (Clarke et al. 2006).This technique allows empty samples to be included in the community analyses.All analyses were performed in R statistical software using packages vegan (Oksanen et al. 2019), ecodist (Goslee and Urban 2007), ggplot2 (Wickham 2009), and grid (R Core Team 2019).

Environmental Conditions
Bloom conditions differed between the study years in three aspects, namely, magnitude, duration, and time between successive hypereutrophic blooms.Heterosigma akashiwo attained a substantially higher maximum density during 2018 (24,996 cells•mL −1 ) compared to pilot study in 2016 (10,395 cells•mL −1 ) (Fig. 3).The number of sampling events between recorded bloom events was greater during the 2016 study, i.e. lasting up to 12 days compared to 3-6 days during 2018 (Fig. 3).Furthermore, the MM site experienced phytoplankton bloom conditions for the entire study period during 2018, but not during the pilot study in 2016.
Bloom presence was more frequent in the MM site and lowest at the lower end of this zone where higher salinity occurs (16-20 salinity) (Fig. 3).Furthermore, while hypereutrophic bloom conditions frequently occurred in the MM and UM sites, they were never recorded in the lower reaches of the mesohaline zone (LM).Accordingly, the density of H. akashiwo in the LM was lower than that observed in the MM and UM sites (Fig. 3).Bloom density was similar in the surface and bottom-waters.A comparison of the diatom:flagellate ratio between years indicated higher ratios during the pilot study 2016 (range 0.04-199.33)than in 2018 (range 0.00-6.45).Additionally, phytoplankton biomass was higher in 2018 (range 24.86-626.78µg chl-a•L −1 ) compared to 2016 (range 3.55-427.42µg chl-a•L −1 ).
The physico-chemical environment was assessed based on four aspects, namely, year, position in water column, site, and bloom presence (Table 1).Physico-chemical conditions were similar during 2016 and 2018.Temperature tended to follow ambient conditions with similar ranges recorded during both years (2016 19.94-26.25 °C, 2018 19.19-25.98°C) (Fig. 4).When comparing surface and bottom waters, temperature ranges were similar and reflected ambient trends (surface 19.81-26.20 ºC, bottom 19.19-26.25 °C) (Fig. 4).Dissolved oxygen was lower at the bottom of the water column, with the maximum recorded concentration of 13.77 mg•L −1 in the surface waters and the lowest of 3.3 mg•L −1 recorded at the bottom.In terms of bloom presence, dissolved oxygen concentration was considerably higher during hypereutrophic blooms (mean concentration 10.33 mg•L −1 ) compared to standard blooms (W = 223.50,P < 0.001) and bloom absence (W = 825.5,P < 0.001).Denser H. akashiwo blooms recorded during    the day positively correlated with higher dissolved oxygen concentrations at night (ρ = 0.343, P < 0.001).Additionally, inverse correlations with dissolved oxygen were observed between salinity (ρ = − 0.264, P = 0.007) and turbidity (ρ = − 0.198; P = 0.045).

Larval Fish Density
A total of 36,711 larvae, belonging to 18 species from 11 families, were collected and identified during the study (Table 2).Gilchristella aestuaria (estuarine round herring) dominated total catch, contributing 96.1% to the total larval density in the mesohaline zone.These were also the only larvae of which all developmental stages were recorded.The most species from a single family came from Gobiidae and collectively contributed 3.4% of the total density observed during the study.Estuarine resident larvae dominated the catch, whereas marine-spawned larvae contributed little to the total density.One freshwater larval fish was caught during the study (family Cyprinidae).Bloom presence did not influence the developmental composition of the larval fish assemblage (Fig. 5).During all three levels of bloom conditions ("none", "standard", and "hypereutrophic"), preflexion stages contributed the most to the assemblage, exceeding 60% under all three bloom conditions.Percentages contribution of preflexion larvae were followed by flexion and then postflexion larvae (Fig. 5).When comparing each developmental stage among bloom presence, no notable differences were found (Table 1).Larval estuarine resident (means: LM = 568 per 100 m 3 , MM = 800 per 100 m 3 , UM = 1311 per 100 m 3 ) and marinespawned density (means: LM = 1 per 100 m 3 , MM = 5 per 100 m 3 , UM = 5 per 100 m 3 ) differed among sites (Fig. 6).Marine-spawned larvae were denser in the UM compared to the MM and LM sites.Estuarine resident larval density was also greater at the UM site compared to the LM and MM (Fig. 6).Estuarine resident larval density did not differ between years (means: 2016 = 792 per 100 m 3 , 2018 = 804 per 100 m 3 ), but marine-spawned density was higher during the pilot study in 2016 (means: 2016 = 7 per 100 m 3 , 2018 = 3 per 100 m 3 ) (Fig. 6).Marine-spawned larvae were also the only group to differ between the surface and bottomwaters, with higher density recorded at the bottom of the water column (means: surface = 2 per 100 m 3 , bottom = 6 per 100 m 3 ) (Fig. 6).
The density of larvae did not differ greatly among bloom severities ("none", "standard", "hypereutrophic").However, the pattern in density showed that the range of density during bloom absence and standard blooms were similar and had higher means (Fig. 6).Smaller ranges, a lower mean, and lower maximum density of estuarine resident larval were recorded during hypereutrophic blooms (range 11-3124 per 100 m 3 , mean 701 per 100 m 3 ) compared to standard blooms (range 4-8444, mean 1146 per 100 m 3 ) and periods of bloom absence (range 2-8715, mean 730 per 100 m 3 ).The density of marine-spawned larvae remained similar across bloom severities.The variation in estuarine resident larval density was best explained by temperature, salinity and dissolved oxygen concentration (AIC = 381.2;deviance explained = 24.3%).Density was negatively related to salinity and dissolved oxygen (Table 3).However, while the GAM indicated that density increased with decreasing oxygen, the density of estuarine resident larvae dropped suddenly at dissolved oxygen concentrations below 5 mg•L −1 (Supplementary Material Figure S2).Marine-spawned larvae were positively related to dissolved oxygen and salinity (AIC = 460.3;deviance explained = 24.2%)(Table 3).As with estuarine residents, there were far fewer marinespawned larvae at dissolved oxygen concentrations below 5 mg•L −1 (Supplementary Material Fig. S3).Marinespawned larval density did not exceed 5 per 100 m 3 at dissolved oxygen < 5 mg•L −1 but ranged from 0 to 50 per 100 m 3 at concentrations > 5 mg•L −1 .

Larval Fish Assemblage Structure
There were no notable differences in diversity among levels of bloom presence (Table 1).While there was overlap in range of diversity, estuarine resident diversity was less varied and had a slightly lower mean (0.19) during hypereutrophic blooms compared to standard blooms and bloom absence (0.22 and 0.21, respectively) (Fig. 5).Similarly, diversity of marine-spawned larval fishes tended to be lower during hypereutrophic blooms (0.08) compared to standard blooms (0.11) and bloom absence (0.15) (Fig. 5).Estuarine resident diversity differed notably among sites (Table 1).The range of estuarine resident diversity was broader at the LM site (0-1.12),whereas the MM and UM sites were similar (0-0.61).In terms of water column position, estuarine resident diversity was greater in surface waters (Fig. 5).In terms of year, estuarine resident diversity was greater in 2018 compared to 2016 while marine-spawned diversity was lower (Table 1).
The estuarine resident larval fish community remained similar in terms of location within the mesohaline zone, water column position, and year (Table 4; Fig. 7).The greatest community difference occurred between the LM and UM sites (R = 0.308) which were best explained by temperature (P = 0.011) and salinity (P = 0.001).There was a noteworthy relationship between the estuarine resident community and H. akashiwo density (P = 0.057).The salinity and H. akashiwo vectors also acted on the community in opposite directions (Fig. 7) and likely acted simultaneously to result in differences among the various sites.The two upper points on the NMDS plot (Fig. 7) were substantially different from all the other samples because they did not contain G. aestuaria and were dominated by Omobranchus woodi and Caffrogobius gilchristi.These samples coincided with the greatest diatom-flagellate ratios (i.e.samples were dominated by diatoms) of 4.79 and 6.45.
Marine-spawned larval fish community remained similar among sites and between years, while only a slight difference was recorded between the surface and bottom (P = 0.023 in Table 4).Salinity (P = 0.012), dissolved oxygen concentration (P = 0.004) and temperature (P = 0.028) contributed to the structure of the marine-spawned larval community (Fig. 8).Empty samples (samples that lacked marine-spawned larvae) occurred across the entire range of temperature, salinity, turbidity, dissolved oxygen, and H. akashiwo densities.Salinity and temperature characterised the community in the surface and bottom waters whereby bottom salinity was higher than surface.In contrast, variation in community structure among sites was related to dissolved oxygen concentration, with the MM and UM sites co-occurring with higher oxygen concentrations.However, there were no distinct differences in the marine-spawned larval community.
The largest difference in the estuarine resident community was among sites (i.e.salinity ranges) (Table 4).Gilchristella aestuaria and Redigobius dewaali were the largest contributors to differences among sites, contributing between 30 and 45% to dissimilarity.Both species were more abundant in the UM sites and decreased in mean density from the upper to the lower end (higher salinity) of the mesohaline zone of the estuary.The estuarine resident larval community did not differ in terms of year or position in the water column.The marine-spawned community only differed between the UM and LM sites.Rhabdosargus holubi, Chelon.dumerili, and Monodactylus falciformis consistently contributed the greatest proportion to the dissimilarity among the various groups (sites, position in water column, and year).There were no C. dumerili recorded during hypereutrophic blooms.
The variability in community structure of estuarine resident and marine-spawned larvae was lower during hypereutrophic blooms (Fig. 9).However, the hypereutrophic blooms was not associated with a clear shift in community structure, whereby R values calculated by ANOSIM was very low (Table 4).As with the differences among sites, water column position, and year, the communities at the various levels of bloom presence were consistently defined by the same dominant species.Differences within the estuarine resident larval fish community were mostly described by G. aestuaria (34-40% dissimilarity) and R. dewaali (27-39%) presence (Table 4).The main contributors to differences among the marine-spawned larval fish community were R. holubi (50-70%) and C. dumerili (10-20%) (Table 4).Larval Fig. 5 Boxplots of estuarine resident and marine-spawned fish larvae density, diversity and species richness compared by bloom presence (NB, no bloom; SB, standard bloom; HB, hypereutrophic bloom); site (LM, lower mesohaline; MM, middle mesohaline; UM, upper mesohaline), depth and year.P-values are adjusted for multiple testing (Bonferroni) ◂ Fig. 6 Developmental stage composition of larval fish assemblages during levels of bloom presence.Note that due to very low density, contributions of leptocephalus larvae and glass eel stage of Anguillidae were not visible on the graph mullet, Chelon dumerilii, and Pseudomyxus capensis were absent during hypereutrophic blooms.

Discussion
Blooms of Heterosigma akashiwo, like many other HABs, altered water quality in the ecologically important mesohaline zone of the Sundays Estuary.In this case, the most apparent effect on the physico-chemical environment were drastic shifts in dissolved oxygen concentration.During the current study, there was a correlation between dissolved oxygen concentration and H. akashiwo density.The high photosynthetic activity synonymous with algal blooms has been demonstrated to substantially increase dissolved oxygen concentrations (Lemley et al. 2018b) which was demonstrated by the high and supersaturated dissolved oxygen concentrations during the current study.These elevated dissolved oxygen concentrations remained present at night because the H. akashiwo bloom had not collapsed (Lemley et al. 2018a).During the current study, phytoplankton biomass was at a near-record high (max.627 µg chl-a•L −1 ) compared to past studies in the Sundays Estuary (e.g.590 µg chla•L −1 ) (Lemley et al. 2017b).Compared to the traditional definition of a bloom (> 20 µg chl-a•L −1 ), such blooms are extremely dense, and it was therefore not surprising that supersaturated oxygen concentrations (up to ~ 13 mg•L −1 at night, ~ 21 mg•L −1 during the day) occurred.Therefore, it appears that HABs of H. akashiwo have the potential to shape larval fish communities by driving DO dynamics.
The overall lack of clear differences in the larval assemblages among various bloom severities may be explained by the estuarine quality paradox (Elliott and Quintino 2007;Borja et al. 2012;Lemley 2015).Estuaries are characterised by a high degree of spatial and temporal variability making these systems highly dynamic (e.g.tides, influx of both marine and fresh water, floods, mouth dynamics).Subsequently, estuarine species are adapted to be tolerant and resilient in the face of major and rapid changes in their environment (Borja et al. 2012).Due to these characteristics, the stress of anthropogenic disturbance on estuarine biota may often go undetected (Elliott and Quintino 2007).
There is little research on the repercussions of oxygen supersaturation on the larval stages of fishes in natural settings.Some studies have concluded that high saturations of oxygen may lead to the formation of bubbles in the guts and mouths of newly hatched larvae, resulting in death (Dannevig and Dannevig 1950;Peterson 1971).Lethargy has been observed in fry populations of Guadalupe Bass Micropterus treculii exposed to hyperoxia (33.3 mg L −1 ) and recommended dissolved oxygen to be kept below 15 mg‧L −1 at 18 ºC (Matthews et al. 2017).Other research on cultured juvenile fishes also recorded gas-bubble disease along with behavioural changes such as reduced swimming (Espmark   There was no consistent pattern in the differences among assemblages by site, depth, year, or bloom conditions.It was expected that density and diversity would be higher at an intermediate level of disturbance (i.e. standard blooms; intermediate disturbance hypothesis) (Connell 1979).When comparing community composition, assemblages overlapped, with the only difference being lower variability during greater bloom presence.Similarly, the pattern in mean density indicated higher abundances during standard blooms as compared to the two extremes (i.e.absence or hypereutrophic).
Dissolved oxygen concentration and salinity were the most consistently important variables during this study as highlighted by GAMs.Salinity zones play an important role in shaping larval fish assemblages in South African estuaries where larval fish density is typically greatest in the mesohaline zone (Strydom et al. 2003(Strydom et al. , 2023;;Strydom 2015).The current study further assessed distribution of larval stages within the mesohaline zone of an estuary.Estuarine resident larval density closely followed that of G. aestuaria during the current study and tended to decrease down the mesohaline zone with increasing salinity.Gilchristella aestuaria typically spawns in the upper reaches of estuaries after which larvae move to the mid-estuary as they develop (Whitfield 1989;Sutherland et al. 2012).
Community structure remained similar to preceding records in the Sundays Estuary (Sutherland et al. 2012).Larval fish species diversity in the mesohaline zone was low and dominated by G. aestuaria (~ 96% of total density).The dominance by G. aestuaria is typical within the mesohaline zone (Harrison and Whitfield 1990;Strydom 2015).While South African estuaries are usually characterised by low diversity and dominance by a few species (Whitfield 1994b), the mean diversities recorded during the current study was below the mean recorded in the mesohaline zone of permanently open South African estuaries (0.51) (Strydom et al. 2003).The second greatest contribution of larvae came from Gobiidae (~ 3%), also commonly found in South African warm-temperate estuaries, although usually more prevalent (Strydom et al. 2003;Strydom 2015).The absence of larval mullet, C. dumerilii and P. capensis, during hypereutrophic blooms is consistent with findings from studies in the Sundays Estuary which reported gill alterations (Bornman et al. 2022a), negative impacts on abundance (Bornman et al. 2022b), and active avoidance (Bornman et al. 2021) by Mugilidae.
Phytoplankton blooms as dense as those recorded during the current study should pose some form of hindrance to larval fishes.The literature on the harmful properties of H. akashiwo highlights three main mechanisms, namely production of toxins, reactive oxygen species, and mucilage (Twiner et al. 2001) which can vary between different strains of the species (Smayda 1998;Higashi et al. 2017;Seoane et al. 2017).It is yet to be determined whether these harmful properties are exhibited by the strain of H. akashiwo that occurs in the Sundays Estuary.Mucilage was observed in conjunction with H. akashiwo blooms during the study period, but its effects remain unknown.Concurrent research on juvenile Mugilidae during this study period, found gill histological changes and lower proportions of secondary lamellae available for gas exchange in fish deliberately exposed to H. akashiwo in situ (Bornman et al. 2022a).Some research has found that H. akashiwo can adhere to the bodies of copepods, effectively impeding swimming and affecting prey ingestion, growth, and survival (Yu et al. 2010;Almeda et al. 2011).In fish, it has been suggested that mucilage may clog gills (Ling and Trick 2010); however, it is not possible to draw such conclusions without assessing health and mortality rates of fish larvae.
Research on other ichthyoplankton in harmful algal blooms has indicated a mix of effects.In the case of a constructed dam, strong harmful algal blooms of Microcystis spp.were linked to low larval fish abundance, whereas reproductive activity remained consistent at low HAB densities (Cataldo et al. 2020).Similarly, the current study found greater peaks in larval density during SB conditions (i.e.lower HAB densities) compared to larval densities recorded during HB conditions.Demirel (2015) reported a positive relationship between larval fish densities and chl-a concentration but that regional-scale trends showed lower larval densities where chl-a concentrations were highest.During Gymnodinium breve blooms outside the inlets of an estuary in North Carolina, lower larval fish densities were recorded although larval recruitment success into the estuary was high later in the season.
The mesohaline zone of the Sundays Estuary still serves as an estuarine nursery to the larval stages of both resident and marine-spawned fish species based on the prevalence of larvae recorded there in this study.However, this highly productive zone now experiences persistent and recurrent HABs (Lemley et al. 2017b) with highly fluctuating oxygen levels.The current study provides some evidence that larval fish assemblages may become less dense and diverse during hypereutrophic blooms in the mesohaline regions of an estuary, possibly driven by the associated dissolved oxygen supersaturation.However, results on this were not very clear and might well be that community level analyses are not good indicators of HAB pressures on larval fishes.A concurrent study on the biochemical body condition of G. aestuaria found that larvae were in poorer condition during hypereutrophic blooms, making them more susceptible to succumb to environmental pressures (Smit et al. 2023).More research is required to fully contextualise the effects of HABs on the early life history stages of fishes.These changes may have major implications for fishes that rely on the Sundays Estuary as a nursery, especially those intricately linked to the estuarine food web relying on G. aestuaria as a prey species.Nutrient pollution from agriculture in regulated, low-inflow estuaries make these systems particularly vulnerable to HAB persistence and subsequent water quality fluctuations, many of which are deleterious for successful grow-out of fish larvae into juveniles.This study demonstrates such a threat and has given insights not only into larval fish dynamics under HAB conditions, but also within a low-inflow estuary, where this information is currently lacking.

Fig. 1
Fig. 1 Geographical location of the Sundays Estuary on the southeast coast of South Africa.Location of sampling area within the mesohaline zone indicated within lines on the estuary (LM, lower mesohaline; MM, middle mesohaline; UM, upper mesohaline)

Fig. 2
Fig. 2 Relationship between H. akashiwo cell density and phytoplankton biomass (chl-a) fitted with a 3-order polynomial trendline (broken line) (equation given on graph)

Fig. 9
Fig.9Non-metric multidimensional scaling plot of estuarine resident (stress = 0.12) and marine-spawned (stress = 0.09) larval fish communities within the mesohaline zone of the Sundays Estuary during spring bloom conditions.Polygons connect all sampling events within

Table 1
Results of the non-parametric analysis of variance (Kruskal-Wallis and Mann-Whitney tests) applied to assess spatial and temporal differences (factors: bloom presence, site, depth, year) in environmental variables, larval fish densities, diversity, and richness (dependent variables)

Table 3
Outputs of statistically selected best models produced through generalised additive modelling (GAM).Covariate statistics include the following: P, approximate significance of smoothed terms; edf, estimated degrees of freedom where edf closer to 1 indicate a linear effect and greater than 1 show a non-linear effect; test statistic: F for estuarine resident models, chisquared for the marine-spawned model; relationship where + indicates an increase in the response variable with an increase in the covariate, and -indicates a decrease in the response variable with an increase in the covariate.Model statistics include the following: AIC, Akaike information criterion (lower value indicates better fit); family, Gaussian for estuarine resident larvae, negative binomial for marine-spawned larvae; link, log; smoothing parameters estimated using marginal likelihood (ML)

Table 4
ANOSIM and SIMPER test results for estuarine resident and marine-spawned communities in the Sundays Estuary over spring periods in 2016 and 2018.R closer to 1 indicates greater dissimilar-