A Snapshot on the Distribution of Coastal Phytoplankton Communities in Five HAB-Affected Bays in Eastern Visayas, Philippines

Introduction

Phytoplankton are known as biological indicators of water and its ecological health status (Wu, 1993; Roelke et al., 2003). They contribute to almost half of the global net primary productivity and are responsible for sustaining the aquatic food webs in support of organisms in the higher trophic levels (Behrenfeld and Randerson, 1998). Phytoplankton respond quickly to different environmental parameters including physical (i.e., light, temperature, rainfall), chemical (inorganic carbon, nutrient, dissolved oxygen), and biological (competition and predation) factors (Maddux and Jones, 1964). These factors have implications in their growth, abundance, availability, and distribution. These understandings on phytoplankton ecology mostly come from studies in temperate areas. Knowing which of these aspects exert a greater influence to the behavior of phytoplankton species in tropical marine waters becomes exceedingly important in understanding the condition of the environment they proliferate, especially in the context of harmful algal blooms. This will also provide knowledge on what to anticipate when seawaters warm due to climate change.

Harmful algal blooms (HABs) are naturally occurring phenomenon that has exacerbated because of anthropogenic disturbances which results to proliferation of microalgae that exhibit a negative impact on aquaculture and aquatic ecosystem or toxic to human health and other aquatic organisms (Van Dolah, 2000; Wells et al., 2020). In Southeast Asia, for example, the Philippines is reported as one of the most impacted regions by HABs (Azanza and Max Taylor, 2001). The first reported case of HAB in the country was in Maqueda Bay, western Samar in 1983 (Maclean, 1989), and since then, the occurrences of these blooms have been increasing spatially and frequently, even diversifying into different HAB types in the recent years (Yñiguez et al., 2020). While many studies about HABs have been reported in the Philippines, most of these are primarily focused in the coastal towns of Bolinao and Anda in the Northern part of the Philippines. Blooms recorded here were caused by different phytoplankton species such as Pyrodinium bahamense (Azanza and Max Taylor, 2001), Alexandrium minutum (Azanza and Benico, 2013; Benico, 2015), Prorocentrum cordatum (Azanza et al., 2005; San Diego-McGlone et al., 2008), Chattonella subsalsa (Lum et al., 2019), Skeletonema sp. (Escobar et al., 2013), Rhizosolenia sp. (Azanza and Benico, 2013), Takayama sp. (Benico et al., 2018, 2019) and many others. Manila Bay, on the other hand, is the socio-economic center in the country draining Metro Manila and many surrounding watersheds also experienced blooms brought about by P. bahamense (Azanza and Max Taylor, 2001) and A. minutum (Benico, 2015). Unfortunately, there has been little published data on the primary hotspots of HAB in the Philippines, which is Eastern Visayas, where Pyrodinium blooms have been experienced frequently.

Region VIII, Eastern Visayas consist of three main islands, Leyte, Samar, and Biliran. This is the site of the first report of toxic Pyrodinium red tides in the Philippines in 1983. And while incidences of Pyrodinium blooms have been increasing in recent years with red tide bans imposed every month (Meniano, 2020; BFAR Region-8, 2021), there has been no phytoplankton ecological studies that have been comprehensively done in the area. This is disappointing as the islands primarily rely on fish, shellfish, and other marine resources as their food and livelihood. P. bahamense incidences of human deaths from mussel consumption have also been reported in recent years. Of the 3,800 PSP cases recorded worldwide during the 1985–2018 period, 2,555 cases of which occurred in the Philippines alone recording a total of 165 deaths (Hallegraeff et al., 2021) of which the majority of these cases originated from western Samar and Biliran Islands. While the Bureau of Fisheries and Aquatic Resources Regional Office 8 (BFAR-8) has been monitoring the HAB occurrences in the region, it has only been limited to HAB species, particularly Pyrodinium bahamense with no physico-chemical data. The successional dynamics of the phytoplankton species in this area has yet to be mapped out. This information can later be used as indicators for prediction of PSP events. This can also serve as seed data toward finding and building physical solutions of HABS in Region VIII. At the same time, we hope to contribute to the understanding of the ecology of not only Pyrodinium, but other phytoplankton species which may or may not cause harmful algal blooms in tropical marine waters.

To obtain the distribution of coast community phytoplankton in Leyte and Samar, a simultaneous monitoring activity was carried out in five major bays in the province of Samar, Leyte, and Biliran in Eastern Visayas from August 2020 to January 2021, which represents the Habagat (Southwest) to Amihan (Northeast) monsoons. The composition, abundance, and distribution of the entire phytoplankton community were determined with focus on microphytoplankton and HABs species. We also correlated the interaction between various physico-chemical parameters to the phytoplankton community and in particular on the toxic species.

Materials and Methods

Sampling Site

The sites chosen in this study were among the bays in Eastern Visayas that experience recent year-round blooms of P. bahamense. In Samar, the three major bays are located within its mussel mariculture zone, namely Cambatutay Bay (11.900 N; 124.786 E), Irong-Irong Bay (11. 818 N; 124.825 E), and Villareal Bay (11.602 N; 124.913 E) were selected. All the bays exhibited a dark-muddy substrate as they received the sediment-rich water from different rivers and other waterway systems. Aside from being muddy, the sediment substrate in Villareal bay was also mixed with granular grains coming from crushed shells. On the other hand, the site in Leyte (Carigara Bay; 11.327 N; 124.693 E) and Biliran (Biliran Bay; 11.556 N; 124.605 E) were beside the coastal town community. There were two stations included in each site resulting in a total of 10 stations (Figure 1).

Figure 1. The map of 10 stations surveyed for phytoplankton composition in Eastern Visayas, Philippines. Sampling sites are designated in different shapes.

Phytoplankton Collection and Analysis

Surface water samples for qualitative analysis were collected monthly in triplicate in each station within 2-m depth using a 20-μm mesh size plankton net with a 30-cm diameter mouth opening. The plankton net was towed vertically and the concentrated 50-mL water was transferred in a translucent bottle and immediately fixed in situ with a 5-mL formalin solution. At the same time another set of surface water samples were collected by a bucket in triplicate in each station within 1-m depth and transferred into a 1-L straight bottle without corrugation to guarantee the settlement of all phytoplankton in the bottom. The samples were immediately fixed in situ using 1% formaldehyde solution and kept in a storage box until the arrival in the laboratory then left undisturbed for at least 24-h. The upper layer of water was carefully sucked out using a plastic capillary tubing and pump and the remaining 250-mL was transferred to a straight container with a relatively small diameter and settled again for another 24-h. The upper layer of water was discarded once again leaving the remaining 50-mL of seawater and the concentrated collected phytoplankton as the final samples. Phytoplankton were quantified using 1-mL aliquot from the final sample that was transferred to a gridded Sedgewick-Rafter counting chamber and allowed to settle for a few seconds before counting at least 300 cells. The densities (cells L^–1) were calculated by multiplying the concentrated volume with the cell counts. Counting was carried out in a 100X total magnification compound microscope in triplicate for every sampling bottle. Micrographs were taken and used in the identification of phytoplankton following the identification guides of Tomas (1997), Larink and Westheide (2006), Omura et al. (2012), and the AlgaeBase database (Guiry and Guiry, 2020). Identification of phytoplankton taxa were limited to diatoms and dinoflagellates.

Physico-Chemical Collection and Analysis

All phytoplankton sampling was accompanied by triplicate measurements of physical (surface water temperature, depth, turbidity, current velocity, and total suspended solids) and chemical (pH, salinity, dissolved oxygen, chlorophyll-a, phosphate, and nitrate) parameters. Temperature, pH, and dissolved oxygen were measured in situ using a calibrated HANNA multiparameter 1-m below the water surface layer. In the same way, salinity was measured using a handheld refractometer. Current velocity was estimated by drift method using a fabricated Holey Sock drogue attached to a 5-m thin rope. The drogue was deployed in the water and the time it traveled within 1-m distance was recorded. Both water depth and turbidity were estimated using a Secchi disc. The depth was obtained by submerging the disc up to the bottom of the water and recording the measurement using the calibration on the rope. Secchi Disk Transparency (SDT; turbidity of the water column) was measured by recording the depth where the disc disappeared from the view (A), and reappeared when hoisted (B), the formula as follows: limit of visibility = (A + B)/2. Data from rainfall was extracted from Philippine Atmospheric, Geophysical, and Astronomical Services Administration (PAG-ASA).

Analyses of water samples for Total Suspended Solids (TSS), chlorophyll-a (chl-a), and dissolved inorganic nutrients such as phosphate (PO₄^–3) and nitrate (NO₃^–) were collected in every sampling period using a 1-L pre-labeled polyethylene (PET) bottles and stored immediately in an icebox. All samples were filtered using 47-mm diameter Whatman^® GF/C glass fiber filters in a motor-driver filtering set-up, while samples for chl-a were done in a low light environment to avoid photodegradation. TSS concentrations were measured following the American Public Health Association [APHA] (1992) method 2540D. Around 250-mL of the collected water was filtered in pre-weighed and pre-dried glass fiber filters. The filters were oven dried at 60°C for 24-h and the total suspended solids was calculated by the difference in initial and final constant weights. The same known volume of water was filtered in a new glass fiber filter for the measurement of chl-a. The filters were soaked in a 10-mL 90% acetone solution at 4°C for 24-h. The samples were then centrifuged for 20 min at 500-g and the absorbance at 750, 663, 645, and 630 nm were measured using a UV-VIS spectrophotometer. The samples were then treated with 100-μL of 0.1N hydrochloric acid (HCl) and the absorbance at 750 and 665 nm was measured once again. Chlorophyll-a concentration was evaluated using Environmental Sciences Section [ESS] (1991) method 150.1. All filtered water samples were then used for nutrient analyses. A volume of 250-mL was transferred in a 150-mL acid-washed beaker and the measurements were provided in full detail by United States Environmental Protection Agency [USEPA] (1978) methods 365.3 (1978) for PO₄⁻³ at 690 nm and PI^® Reef Master Test Kit for NO₃^– at 540 nm using the spectrophotometer. The collection of water samples was done in triplicate in each station, and at least three replicate measurements were carried out in each analysis.

Statistics

Data for phytoplankton abundance and physico-chemical parameters were averaged monthly in each site to acquire mean semi-annual data. The mean abundance data were utilized to obtain biodiversity indices such as Shannon-Wiener diversity (H′), Margalef’s richness index (D), and Pielou evenness (J′). using the Paleontological Statistics (PAST) software program version 4.13 (Hammer et al., 2001). The Principal Component Analysis (PCA) was carried out to evaluate the phytoplankton community and physico-chemical parameters at different bays and the extent to which the measured parameters could explain the distribution of the phytoplankton species. Before inputting into PCA, the mean data were centered and/or scaled to normalize data that were associated from different magnitudes followed by standardization for the physico-chemical parameters or pre-transformation using Helinger for the phytoplankton abundance. The PCA was carried out using the prcomp function. The scores for PCA 1 and 2 were extracted and used to visualize the PCA biplot of physico-chemical parameters using the ggbiplot2 (Wickham, 2016), phytoplankton taxa using factoextra (Kassambara and Mundt, 2017), and dominant phytoplankton species against the environmental parameters using Factoshiny package (Vaissie et al., 2021). To assess the correlation between the environmental drivers and dominant phytoplankton genera in each bay, the Pearson’s correlation coefficient was performed using the corr function. We then performed a correlation test to know whether the correlation coefficient between variables is significantly different from 0 using the rcorr function. A correlogram was built using corrplot (Wei and Simko, 2021) to visualize correlation between all possible pairs present in the dataset, with a clear distinction for correlations that are different from 0. All multivariate analyses were performed in R.

Results

Physico-Chemical Profiles in Eastern Visayas Sites

Monthly time-series for the physical and chemical environmental data recorded during the sampling period for the five HAB-affected bays in Eastern Visayas are presented in Figure 2. On the average, the surface water temperature in all the bays were collectively high at 29–32°C. The highest peak in temperature occurred in September in all the bays with the highest value recorded in Cambatutay Bay at 32°C, however, no obvious pattern is apparent in the temperature plot. The same is true in the pH values that were relatively alkaline (7∼8) and did not very much in between the five bays. Interestingly, the three bays in Samar seem to follow the same pattern compared to the bays in Leyte and Biliran in terms of turbidity, salinity, current velocity, dissolved oxygen, and TSS. The average depth in Cambatutay, Villareal, and Irong-Irong Bays were all shallow at 4-m, high in turbidity and current velocity, but low in salinity (Cambatutay = 30 ∼ 34 ppt; Irong-Irong = 31 ∼ 35 ppt; Villareal = 27 ∼ 34 ppt), dissolved oxygen (Cambatutay and Irong-Irong = 5 ∼ 7 ppm; Villareal = 5 ∼ 6 ppm), and TSS (Cambatutay, Irong-Irong, and Cambatutay = 0.3 ∼ 0.1 mg/L) compared to Carigara (salinity = 30–36 ppt; dissolved oxygen = 8 ∼ 9 ppm; TSS = 0.02 ∼ 0.2 mg/L) and Biliran (salinity = 30 ∼ 36 ppt; dissolved oxygen = 8 ∼ 9 ppm; TSS = 0.03 ∼ 0.02 mg/L) Bays. According to PAG-ASA, the highest rainfall during the sampling period starts in October in Leyte while in Samar, the increase in rainfall starts in December. Yet on the average, the rainfall concentration remains relatively high (Leyte = 162 ∼ 692 mm; Samar = 251 ∼ 556 mm) in Eastern Visayas. The alternate level on the concentration of chl-a was observed monthly for most of the bays, however, the average chl-a level in Carigara and Biliran Bays were a little high compared to Irong-Irong, Cambatutay, and Villareal bays. The PO₄^–3 concentration was relatively the same in all of the bays at the mean of 0.01 μM. The nitrate concentrations in the three bays in Samar were slightly higher, at 0.06 μM than Carigara and Biliran Bays at 0.04 μM. Although starting in November, the NO₃^– concentration in Irong-Irong, Cambatutay, and Villareal Bays decreased to 0.04 μM.

Figure 2. Average measurement of physico-chemical parameters in the five HAB-affected bays in Eastern Visayas from August 2020-January 2021, including the inorganic nutrients (nitrates and phosphates). The number of parameters measured varied depending on the availability of the equipment during sampling and weather or sea condition.

Principal Component Analysis (PCA) using the physico-chemical parameters of the five bays used in this study separated the samples into two distinct clusters (Figure 3A), which explained 44.2% of total variance by the first two principal components (PCs): PC1 = 26.1%; eigenvalue = 3.13 and PC2 = 18.1%; eigenvalue = 2.17. The bays of Cambatutay, Irong-Irong, and Villareal clustered together and were separated against the bays of Carigara and Biliran. Looking at the axis, the bays for Biliran and Carigara were closely correlated on salinity, pH, dissolved oxygen, depth, chl-a, and PO₄^–3 while Cambatutay Bay was closely correlated with rainfall and current velocity.

Figure 3. Biplot of Principal Component Analysis (PCA). (A) PCA biplot for physico-chemical parameters of the five bays against different environmental gradients. The axis scores were shown for PCA 1 and PCA 2. There are five groups included during the PCA, represented by the five different bays designated by a specific color. Physico-chemical parameters legend: Sal (salinity), Turb (turbidity), Temp (temperature), Rfall (rainfall), CV (current velocity), TSS (total suspended solids), chla (chlorophyll-a), PO₄⁻³ (phosphate), NO₃^– (nitrate), and DO (dissolved oxygen). (B) PCA biplot of the phytoplankton abundance in five sampling sites. The blue arrows represent the phytoplankton genera. Phytoplankton legend: Ske (Skeletonema), Tha (Thalassionema), Pseu (Pseudo-nitzschia), Cha (Chaetoceros), Cos (Coscinodiscus), Bact (Bacteriastrum), Bac (Bacillaria), Dit (Ditylum), Cer (Ceratium), Noc (Noctiluca), Pyr (Pyrodinium), Det (Detonula), Rhi (Rhizosolenia), Cyl (Cylindrotheca), Ple (Pleurosigma), Dip (Diploneis), Nitz (Nitzschia), Coc (Cocconeis), and Pro (Prorocentrum).

Phytoplankton Community in the Five Bays

A total of 43 phytoplankton genera were identified from the five bays in Eastern Visayas, mainly dominated by the taxa belonging to Bacillariophyceae (33 diatom genera), and nine Dinophyceae (dinoflagellates). The Samar bays had higher monthly total average phytoplankton densities at 55,000–73,000 cells/L compared with Leyte bays at 26,000–31,000 cells/L (Figure 4). For the whole sampling period, Irong-Irong Bay had the highest mean density of phytoplankton species (38,000 cells L^–1, H′ = 1.81, J′ = 0.32, D = 1.83) then followed by Villareal (37,000 cells L^–1, H′ = 0.155, J′ = 0.26, D = 1.83), Cambatutay (32,000 cells L^–1, H′ = 1.91, J′ = 1.78, D = 1.96), Carigara (25,000 cells L^–1, H′ = 1.67, J′ = 0.28, D = 1.92), and Biliran Bays (22,000 cells L^–1, H′ = 1.58, J′ = 0.27, D = 1.83), respectively. No clear spatial pattern can be discerned as to the evenness values of the five bays that remained low throughout the sampling period (0.2 to.4), yet the species diversity of these bays was fairly similar ranging from 1.6 to 1.9. Likewise, the Margalef’s richness index of the five bays was similar with a range of 1.8 to 2.0.

Figure 4. Mean relative abundances (%) of the dominant phytoplankton taxonomic groups found in each sampling site in Eastern Visayas, Philippines from August 2020-January 2021 and the corresponding total average phytoplankton densities in 10³ cells L^–1.

The PCA showed that the first two principal components accounted for more than 71% of the variations (Dim1 = 42.1% with eigenvalue = 7.99 and Dim2 = 28.8% with eigenvalue = 5.47) in the species composition, with obvious pattern of species division between sites (Figure 3B). Top dominant genera, Skeletonema, Pseudo-nitzschia, and Noctiluca scintillans all exhibited positive linear correlation in Villareal Bay while Pyrodinium bahamense displayed positive correlation with Cambatutay and Villareal Bays. Most species with positive correlation in Carigara and Biliran Bays were non-HAB diatoms such as Rhizosolenia, Bacteriastrum, and Chaetoceros.

Based from the phytoplankton community composition, the relative abundances of diatoms, which accounted for up to 99% in Biliran Bay, 95% in Villarreal Bay, 94% in Carigara Bay, and 89% in Irong-Irong and Cambatutay Bays, remained high throughout the sampling period across all the five bays in contrast to dinoflagellates that was collectively < 11% in the population. Diatoms in the genus Skeletonema, Bacteriastrum, Chaetoceros, Pseudo-nitzschia, and Thalassionema were common and prevalent in all the bays during the sampling period (Figure 4). For dinoflagellates, genera which are HAB vectors such as genus Ceratium, Pyrodinium bahamense, and Noctiluca scintillans were common and prevalent in the three bays in Samar. In fact, we observed a spike in the density of Pyrodinium within September to November while the density on N. scintillans starts to increase in December. On the other hand, we observed the presence of other potential HAB vectors such as Dinophysis, Alexandrium, Nitzschia, Prorocentrum, and Gymnodinium, but all were collectively below 100 cells L^–1. Between the dominant genera, the most prevalent group belonged to Skeletonema with relative mean density of 10,000 cells L^–1 but abundances of Pseudo-nitzschia species and N. scintillans were relatively high too, especially in the three bays in Samar, with mean densities of 4,000 and 600 cells L^–1, respectively.

Diatom Succession in the Five Bays

For Irong-Irong Bay in the Samar Seas, the diatoms Skeletonema and Pseudo-nitzschia dominated the phytoplankton community in August 2020, the start of the sampling period. Skeletonema then increased further and dominated in September, to be briefly replaced by Bacteriastrum the month after and Rhizosolenia by the next month, November. The profile in December changed with Skeletonema, Bacteriastrum, and a growing number of Chaetoceros as the most prevalent. Pseudo-nitzschia continued to decrease while Chaetoceros increased further by the end of the sampling period, in January 2021. In Cambatutay Bay, Skeletonema dominated the phytoplankton community except in December when Ditylum, Bacteriastrum, and Chaetoceros briefly dominated. Chaetoceros continued to increase in number in January. Villareal Bay was dominated by Pseudo-nitzschia sp. during August but they gradually decreased in the months after. Unlike Pseudo-nitzschia sp., the density of Chaetoceros increased monthly, with the highest increase in December and January. The abundance profile of Bacteriastrum was heterogeneous but they are mostly abundant in October and December. The density of Skeletonema spiked in September to November to abruptly decrease to low in December then increased by the next month, January.

In the Leyte side (Figure 4), more diatoms dominated in the two bays examined. In Carigara, Chaetoceros dominated the months of August and September. This was briefly replaced by Thalassionema. Then Skeletonema dominated the rest of the months. Chaetoceros also dominated the phytoplankton community of Biliran Bay except for the month of September when this was dominated by Bacteriastum.

Correlation of Harmful Algal Blooms Species and Environmental Drivers

Microscopic images of the HAB vectors identified in the sampling sites are shown in Figure 5. The relationship between the dominant phytoplankton groups with the physico-chemical parameters and potential HAB vectors in the five bays were further subjected to a correlation analysis through Pearson’s correlation coefficient (Figure 6). In Villareal Bay, PO₄^–3 (r = 0.47; p ≥ 0.999), ${\text{NO}}_3^{-}$ (r = 0.48; p ≥ 0.999), and rainfall (r = 0.62; p ≥ 0.999) displayed positive correlation with P. bahamense. The same can be observed in Irong-Irong (${\text{PO}}_4^{-3}$: r = 0.61; p ≥ 0.999; ${\text{NO}}_3^{-}$: r = 0.47; p ≥ 0.999; rainfall: r = 0.13; p ≥ 0.999) and Cambatutay (${\text{PO}}_4^{-3}$: r = 0.4; p ≥ 0.999; ${\text{NO}}_3^{-}$: r = 0.52; p ≥ 0.999; rainfall: r = 0.47; p ≥ 0.999) Bays. Moreover, the relationship of temperature against dominant genera was highly heterogeneous with no clear spatial pattern, however, it displayed a positive correlation with N. scintillans in Villareal bay (r = 0.55; p ≥ 0.999). Pseudo-nitzschia which was partly prevalent in the three bays in Samar displayed positive correlation with dissolved oxygen (Villareal: r = 0.78, p ≥ 0.999; Irong-Irong: r = 0.33, p ≥ 0.999; Cambatutay: r = 0.95, p ≥ 0.999) but exhibited a negative correlation in Carigara (r = 0.07, p ≥ 0.999) and Biliran (r = −0.29, p ≥ 0.999) Bays. Between species interaction, P. bahamense displayed positive correlation with N. scintillans in Irong-Irong (r = −0.59, p ≥ 0.999) and Cambatutay (r = −0.68, p ≥ 0.999) Bays but exhibited a negative correlation in Villareal Bay. Pseudo-nitzschia that was prevalent in all the bays displayed a heterogenous relationship against other genera, except in Villareal (r = −0.50; p ≥ 0.999) and Cambatutay (r = −0.46; p ≥ 0.999) Bays where they exhibited a positive correlation with P. bahamense. Other diatom genera displayed a highly heterogeneous relationship between potential HAB vectors and non-HAB phytoplankton. Other correlations observed between the physico-chemical parameters, dominant phytoplankton taxa, and potential HAB-vectors were not significant.

Figure 5. Selected micrographs of the phytoplankton observed in the five bays of E. Visayas, 2020–2021: (A) Pseudo-nitzschia sp. (B) Prorocentrum sigmoides, (C) Ceratium furca (Tripos furca), (D) Pyrodinium bahamense, (E) Gymnodinium catenatum, (F) Noctiluca scintillans, (G) Dinophysis caudata, (H) Rhizosolenia fallax, (I) Ceratium trichoceros (Tripos trichoceros), (J) Ceratium inflatum (Tripos inflatus), and (K) Thalassionema nitzschioides. Scale Bars: (A–D,F–H) = 10 μm; (E, I–K) = 100 μm.

Figure 6. Correlogram of the dominant phytoplankton groups, physico-chemical parameters, and dominant HAB-vectors in (A) Cambatutay, (B) Irong-Irong, (C) Villareal, (D) Carigara, and (E) Biliran Bays. The correlogram shows correlation coefficients for all pairs of variables. Positive correlations are displayed in blue and negative correlations in red. The intensity of the color is proportional to the correlation coefficient so the stronger the correlation (i.e., the closer to –1 or 1), the darker the boxes. The color legend on the right-hand side of the correlogram shows the correlation coefficients and the corresponding colors. The physico-chemical parameter legend: DO (dissolved oxygen), CV (current velocity), and TSS (total suspended solids).

Discussion

Spatial difference in terms of the physico-chemical profiles was observed between the five bays in Eastern Visayas. The bays located in Samar displayed a different profile than the bays in Leyte and Biliran. This is not surprising as Cambatutay, Irong-Irong, and Villareal Bays are sites of mussel mariculture zone near the river mouth in Samar unlike Carigara and Biliran Bays which are beside coastal communities. It is most likely that the differences observed in these bays will play a significant role in the abundance, distribution, and diversity of phytoplankton communities in these areas (Reynolds, 1984; Klausmeier and Litchman, 2001; Clegg et al., 2004).

The profiles we initially recorded in the bays in Samar were comparable to the study reported by Cebu and Orale (2017). The low concentration of dissolved oxygen and salinity can be attenuated due to heavy rainfalls and the possible volume of freshwater input from the rivers (Damotharan et al., 2010; Fatema et al., 2014; Cebu and Orale, 2017). Moreover, high levels of turbidity may have prevented the light from penetrating in the water column and may, in part, affect the chlorophyll-a concentrations. Since the three bays in Samar are located near the mussel farming area, the decline in concentration of dissolved oxygen may, in part, be associated with how mussels control the water chemistry within its habitat, maintaining an acidic (Martinez Rodriguez et al., 2015), anoxic (Nicklisch et al., 2016), and ionically sparse (Yu et al., 2011; Miller et al., 2015) environment to be able to keep producing byssal threads, a proteinaceous fiber that allows the mussel to adhere to surface underwater thereby preventing dislodgement (Bell and Gosline, 1996; Carrington et al., 2009). Furthermore, mussel farms in the area are extensive, and the farmers practiced the “staking” or “wigwam” method that uses bamboo that is stuck in the sediments. Most of the time these are left in the same area after harvest and thus contribute to the organic matter in the sea. Meanwhile, the concentrations of dissolved inorganic nutrients measured in all of the bays are generally low, an indication of impending or post bloom phase as several red tide outbreaks were announced throughout the duration of the sampling period.

Diatoms are a major component of the biological community, serving as one of the primary oxygen sources in aquatic environments. In Eastern Visayas, we observed the prevalence of diatom associated groups in its phytoplankton community. Among the diatom taxa, the genus Skeletonema was the most prevalent. This result was also observed by Azanza and Miranda (2001) in Manila Bay. The same group has caused fish kills in Bolinao, Pangasinan (Escobar et al., 2013). Other studies elsewhere, such as in the coastal and estuarine areas in China, found that Skeletonema costatum to be a cosmopolitan phytoplankton in the area (Huo and Shu, 2005). In the three bays in Samar, Skeletonema was positively correlated with temperature, salinity, and nitrates, and thus were abundant when these parameters were higher in these bays. This group is also the most prevalent in the waters of Biliran and Carigara Bays and was correlated to nitrate concentration. It has been shown both in the laboratory and in the field that Skeletonema blooms in higher temperatures (Shikata et al., 2008). Moreover, S. costatum has also been reported to uptake large concentrations of nitrate that results in chain elongation (Takabayashi et al., 2006). These results indicate that the waters in the five bays provide a suitable condition, at least with parameters such as temperature, salinity, and nutrients, for Skeletonema to grow in abundance (Zohdi and Abbaspour, 2019). Meanwhile, other diatom associated groups such as Chaetoceros, Bacteriastrum, Pseudo-nitzschia, and among others which are dominant in other parts of the Philippines (Yap-Dejeto et al., 2008; Yap-Dejeto et al., 2013; Yñiguez et al., 2020), were also abundant in all of the bays herein.

Aside from diatom groups, dinoflagellates such as P. bahamense, N. scintillans, and species of Ceratium co-dominated especially in the waters in Samar Sea. N. scintillans and species from genus Ceratium have been reported to cause fish kill (Mijares et al., 1985; Pullin et al., 1993; Orellana-Cepeda et al., 2004; Baek et al., 2008) while P. bahamense are STXs-producing causing Paralytic Shellfish Poisoning. The co-dominance of these dinoflagellate groups along with other diatoms in the three bays in Samar can be associated with the interplay of the environmental factors present in the area. The three bays both exhibited high rainfall and nitrate concentrations. Most of these species also displayed positive correlation with these parameters in these bays. Since nitrate is the primary limiting nutrient for marine phytoplankton (Thomas, 1966, 1969) and heavy rainfall and winds are among other features contributing to HAB development (Mallin and Corbett, 2006), then rainfall-induced nutrient fluctuations is one of the key elements contributing in the co-dominance of HAB and non-HAB phytoplankton in the Samar Sea.

Moreover, the assemblages of potential HAB vectors such as Pyrodinium, Noctiluca, Ceratium, and Pseudo-nitzschia observed in the bays in Samar, of which the sites are mussel farming areas, were also observed by Albelda et al. (2019) in the areas near and within the fish farming site in Bolinao, Pangasinan. There is some evidence indicating that sites for mariculture or wild harvest of shellfish are hotspots of HAB events in the Philippines (Villanoy et al., 2006; Yñiguez et al., 2018). For example, P. bahamense is the causative agent of the red-tide outbreak in Samar. In fact, within September to November, a noticeable increase in the relative density of P. bahamense was observed in the three bays in Samar which also coincided with several red tide bans announced by BFAR 8 on the same duration in the province. In addition, mussels are naturally filter feeders and can accumulate high levels of STXs-producing phytoplankton. It could be that the feces and digestive tract of the mussels from the infected area can be loaded with viable P. bahamense cells and its resting cyst which, when transferred, can act as vector or inoculum to another mariculture site or to the neighboring waters (Bricelj et al., 1991; Hégaret et al., 2008). However, these hypotheses need to be confirmed.

Interestingly, we found a slight (statistically insignificant) positive correlation between P. bahamense and the unarmored dinoflagellate, N. scintillans, in Cambatutay, Irong-Irong, and Villareal Bays. The green N. scintillans was the cause of a fish kill event in Manila Bay in 1987. In 2010, Azanza et al. (2010) reported a prey-predator relationship between P. bahamense and N. scintillans using a feeding experiment in a laboratory setting. In 2017, Folio and Yap-Dejeto (2022) observed a decline of P. bahamense density with increasing N. noctiluca cells in Irong-Irong Bay. It is most likely that the abundance of N. scintillans in the Samar Sea can be associated with the presence of P. bahamense in its surrounding waters, however, enough evidence and a thorough study will be needed to unravel the complexity of this relationship.

In September, we detected an abrupt increase in the density of Bacteriastrum in Carigara and Biliran Bays while this observation occurred in October in Irong-Irong and Villareal Bays. The shifts in Irong-Irong and Villareal Bays must have been caused by elevated nitrates during that month in these bays. This is supported by positive correlations of nitrates with Bacteriastrum in all Samar bays. In Villareal Bay, elevation of nitrates supported other diatoms to bloom, and thus the bloom of Bacteriastrum is shared with Skeletonema. We can only surmise a similar occurrence in the bloom of Bacteriastrum in Carigara and Biliran Bays since nutrient data was not acquired during that time. Additionally, it may have been caused by peaks of temperature in September in both bays and a peak of current velocity in Carigara Bay. Bosak et al. (2016) also observed an abundance of Bacteriastrum at higher temperatures. Higher current velocity could have caused turbulence and water mixing increasing the availability of nutrients to chain forming Bacteriastrum causing these to multiply. Carigara Bay’s Bacteriastrum bloom is only secondary to Chaetoceros bloom which is another chain forming diatom whose nutrient uptake is supported by turbulence (Dell’Aquila et al., 2017).

Chaetoceros dominated the bays of Irong-Irong and partly in Cambatutay Bays during the end of the sampling period (January 2021), while these were abundant during the beginning of the sampling period, in August and September, 2020 in Carigara Bay. In Villareal Bay, these dominated in both the end and the beginning of the sampling period and persisted during the whole of the sampling period. In the coastal area of northeastern Adriatic, slight increase of nutrients, mainly phosphate, influenced the Chaetoceros blooms (Bosak et al., 2016). In this study, Chaetoceros positively correlated with both phosphate and nitrate in the three bays in Samar. But the peaks of nutrients do not match the months when Chaetoceros increased. Current velocity which could affect turbulence is also correlated to Chaetoceros density in the bays where Chaetoceros bloomed. Thus, as Dell’Aquila et al. (2017) suggested, we hypothesize that the increasing densities of Chaetoceros could be an interplay of increased nutrients and turbulence during these months in these bays. The persistence of Chaetoceros in Villareal Bay for the whole sampling period indicates the existence of Chaetoceros resting spores (Ishii et al., 2011) in the area. The absence of these spores might also be the reason why this never bloomed in Biliran. This is just a hypothesis and will need to be verified by further research that will include searching for resting spores in sediments.

The dominance of the bloom forming diazotroph diatom (Villareal, 1992), Rhizosolenia, in Irong-Irong Bay in November must have caused an increase of dissolved oxygen as the data for these two components positively correlated. Other parameters showed only slight correlations and the bloom of this phytoplankton was not observed in any of the other four bays or in any other months during the study period. For example, Yoshimatsu et al. (2020) showed that maximum growth of Rhizosolenia setigera occurs in relatively warmer temperatures. Carstensen et al. (2004) also recorded summer blooms of Rhizosolenia. The temperatures in all five bays were relatively warm. Moreover, we also observed that there was a drop in salinity in Irong-Irong Bay during the time of the bloom. We hypothesize that the decrease in salinity in the month of November in this bay could have caused Rhizosolenia to dominate. Species of Rhizosolenia were noted to be stenohaline or even preferred lower salinities (Gnatiades and Smayda, 1970; Rijstenbil, 1987; Yoshimatsu et al., 2020). This then gave an advantage of Rhizosolenia to bloom over the rest of the other phytoplankton.

Thalassionema, on the other hand, must have provided for the high chlorophyll-a in October in Carigara Bay as both factors peaked during this month. Similar to Rhizosolenia, Thalassionema dominance was only observed once, and only in one bay, Carigara Bay. In our study, Thalassionema negatively correlated with current velocity. But in the study in India, Thalassionema bloomed together with other phytoplankton under turbulent conditions during the Southwest Monsoon (Retnamma et al., 2020). Seasonal change brought about by the Southwest Monsoon could also be the reason for the transient succession of Thalassionema in this study as October is the end month of the Southwest Monsoon in the Philippines. The effect of the monsoon conditions must have influenced the phytoplankton community in Carigara Bay during that time. Canini et al. (2013) also observed Thalassionema bloom during the Southwest Monsoon in Panguil Bay, Philippines.

In December, we detected a sudden increase of density between the three genera, Skeletonema, Bacteriastrum, and Chaetoceros in the Samar Sea. The prevalence of these species was accompanied by the increasing rainfall that started in the same month and subsequent decrease in the nitrate concentration by twofold. The dominance of these species was also evidenced by the increase of chlorophyll-a concentration in the three bays in Samar during December. These results are consistent with our observation from previous months showing that physical disturbances of the water column brought about by heavy rainfall indeed provides a favorable condition for phytoplankton species to grow, and sometimes, blooms. It seems that environmental factors, particularly rainfall, nitrate, and chlorophyll-a are likely associated with the structuring of these phytoplankton species in the Samar Sea. Interestingly, we observed a succession between Chaetoceros and Pseudo-nitzschia in the waters of Villareal, Irong-Irong, and Cambatutay Bays. In the three bays, Pseudo-nitzschia displayed a decreasing density in contrast with the increasing density of Chaetoceros. Unlike Chaetoceros, the Pseudo-nitzschia sp. in the three bays were correlated with decreasing temperature. This result indicates that the successional pattern discerned between the two diatoms might be governed by temperature as a factor driving the structuring of these species. This result is in line with previous observations showing temperature as one of the fundamental factors that influence the structuring of microalgal communities (Eppley, 1972; Karentz and Smayda, 1984; Bouman et al., 2003). Moreover, this data somewhat agrees with studies in the temperate regions that showed Pseudo-nitzschia to be temperature dependent (Lundholm et al., 1997).

Given these insights, we can say that the linkages of the phytoplankton groups displayed here provides a preview on the succession of phytoplankton communities in the waters of the five major HAB-affected bays in Eastern Visayas. The shifts of the phytoplankton composition observed in this study were governed by a myriad of complex interactions of many environmental factors (Barbosa et al., 2010; Pulina et al., 2012). However, there are other strategies that could also influence the turnover in the composition and abundance of phytoplankton such as stratification, grazing, predation, or allelopathy (Čalić et al., 2013). It is equally important to understand the impact of these biological drivers on the phytoplankton community structure and how this component relates to the HAB problem in Eastern Visayas. This emphasizes the need for a more comprehensive study involving many factors to better understand the dynamics in the change of phytoplankton community and determine which of these factors have influenced such changes. Furthermore, these results make the waters in Eastern Visayas a natural mesocosm to test fundamental questions related to tropical marine phytoplankton ecology and phytotoxins. Insights and answers will greatly benefit not only the ongoing problem of harmful algal bloom in the region where potential policy can be crafted and fine-tuned in the need of the area; but it provides a peak of what is to come when waters become warmer and more acidic due to climate change.

Conclusion

The present investigation explores the phytoplankton community structure and the physico-chemical profiles of five HAB-affected bays in Eastern Visayas from August 2020 to January 2021. Based on our initial results, the physico-chemical profiles of Irong-Irong, Cambatutay, and Villareal bays were different compared to the bays in Carigara and Biliran. The distinctness observed between these bays provides insights on the structure of the phytoplankton community present in these areas, especially on species that cause HAB. Trends of physico-chemical parameters as well as phytoplankton succession were documented and discussed. All five bays were diatom-dominated and succession between diatoms varied in the four bays. The diatoms that dominated in the bays were Skeletonema, Pseudo-nitzschia, Bacteriastrum. Chaetoceros, Rhizosolenia and Thalassionema. The following are the HAB causing species identified in the study; Pseudo-nitzschia spp., Nitzschia spp., Alexandrium sp., Pyrodinium bahamense, Gymnodinium sp., Dinophysis caudata, Dinophysis miles, Prorocentrum lima, Skeletonema spp., Noctiluca scintillans, Prorocentrum sigmoides, Ceratium furca, Ceratium fusus, Ceratium inflatum, and Prorocentrum micans. Overall, we see a snapshot on the phytoplankton ecology and HAB dynamics in the five HAB-affected bays in Eastern Visayas. These results will serve as a baseline information for future studies that aims to understand further the problem related to harmful algal blooms in the region.

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.

Author Contributions

LGYD, MLLA, and MBA conceptualized the study. SFR, MLSS, MLLA, JAO, and EGA conducted the data collection. SFR and LGYD conducted analysis and interpretation of the results. SFR, MLSS, and LGYD wrote the manuscript. All authors drafted, approved and agreed to be accountable for all aspects of the work in this article.

Funding

This work was supported by the DOST-GIA project of Department of Science and Technology—Philippine Council for Agriculture, Aquatic, and Natural Resources Research and Development (DOST-PCAARRD) entitled “WATCH HAB R8: Warning And Technology to Combat Hazards of Harmful Algal Blooms in Region 8” to LGYD.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to acknowledge the assistance of Alessandra M. Malacora in the phytoplankton identification and Clances Joanne Urmeneta in phytoplankton microscopy. The authors would also like to thank the Bureau of Fisheries and Aquatic Resources (BFAR) Region 8 for their assistance from the formulation of the study to the conduct and field experimentation. The authors also grateful to the reviewers who were instrumental in tremendously improving the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.730518/full#supplementary-material

References