Dietary Changes in the Ark Clam (Anadara kagoshimensis) Is Associated with Phytoplankton Community Patterns in a Temperate Coastal Embayment

Abstract

The monthly phytoplankton communities and dietary items of the filter-feeding ark clam (Anadara kagoshimensis) in cultivation were examined in a shallow temperate coastal embayment of Yeoja Bay in Korea, to identify dietary changes in clams associated with phytoplankton community patterns. A self-organizing map (SOM) algorithm was applied to shape the community structures of phytoplankton. Clam δ¹³C and δ¹⁵N values were determined monthly and compared with those of phytoplankton, microphytobenthos, suspended particulate organic matter (SPOM), sedimentary organic matter (sedimentary OM), and Phragmites australis. Our SOM clustered monthly phytoplankton communities, revealing a seasonal shift in the dominance of large-sized diatoms (sporadically together with dinoflagellates), which were detected almost year-round, to small-sized taxa (chlorophytes, prasinophytes, and prymnesiophytes), which were observed in May–June. The δ¹³C and δ¹⁵N measurements revealed that pelagic and benthic diatoms serve as the main contributors to the clam diets. A reduction in their dietary contribution accompanied a considerable increment in the contribution of Phragmites detritus in association with the dominance of small-sized phytoplankton during the late spring. Our results suggest that the dominance of small-sized phytoplankton during the critical spring period of the clam life cycle may decrease the availability of preferred items (i.e., size-related food quality) and lead to dietary changes in the clams in relation to climate forcing in this warming sea.

1. Introduction

Shallow coastal embayments are characterized by complex hydrodynamic regimes, sediment dynamics, and biochemical processes [1,2,3]. The embayment ecosystem serves as a convergence of organic matter and nutrients driven by groundwater and surface water discharge, tidal exchange, and coastal upwelling [4,5,6]. The sources of biogenic compounds originating from the embayment ecosystem itself are also highly diverse, including microbial mats, microalgae, and macrophytes [7,8]. The abundance and diversity of the primary sources of organic matter contribute to the variation in the food-web structure [9]. The composition and abundance of primary producers keep changing over time, thus leading to the development of trophic interactions and consequently establishing the structural and dynamic properties of the food web in the embayment ecosystem. The resource–consumer relationships can be characterized by the abundance/biomass of organisms and the resource contribution to the diets of consumers [9,10,11]. The patterns of prey utilization by consumers in response to the varying availability of resources are an ecosystem attribute that determines the strength and length of trophic links between the major food-web components and the embayment ecosystem processes.

The coastal embayment ecosystem is being widely used for bivalve cultivation for human consumption [1,8,9,11]. Those bivalves employ a suspension-feeding mode that filters particles in a water column, and their growth and production are largely dependent on phytoplankton and other trophic subsidies (e.g., microphytobenthos and marsh detritus) [12,13,14]. Because phytoplankton has a high sensitivity to environmental changes and constitutes the base of aquatic food chains [15,16], changes in its communities may affect trophic interactions, which in turn results in modifications in the structure and processes of the entire food web of the embayment ecosystem [9,17,18]. As indicated by the effective pre-ingestive selection of food adopted by many bivalve species [19,20], seasonal and spatial patterns in phytoplankton community structure may facilitate a clear distinction in the resource–consumer interactions across an environmental gradient [21,22,23]. Conversely, a top-down regulation mediated by the bivalves may govern the community structure of prey based on their preferential selectivity according to size and nutritional value [24,25]. Therefore, patterning in phytoplankton community composition will provide a basic feature to delineate the resource–consumer interactions and ecosystem processes associated with the wide environmental conditions that are found specifically in temperate coastal embayment systems and sustain high biological diversity and production.

The classification of community structures based on the composition and abundance of occurring assemblages according to the type of water body at different space and time points, e.g., [21,23] requires the evaluation of the resource contribution to consumers’ diets needs, to identify the role of actual prey items as individual food-chain supporters. The identification of the diets of animals to delineate the resource–consumer relationships cannot be readily performed in natural aquatic systems because of methodological limitations according to quick digestion after feeding, difficulties in identification of too many taxonomic groups, and often empty stomachs [26]. The stable isotopes of carbon and nitrogen (δ¹³C and δ¹⁵N) have been widely applied to quantify the relative contribution of various organic matters to consumer diets in different nodes of food web networks [27,28]. Dietary carbon can be traced in an animal based on its tissue δ¹³C values, because the distinct values between primary producers are conserved during trophic transfers, with little or no trophic enrichment in ¹³C (≤1‰) [29]. In turn, δ¹⁵N values are used to estimate the trophic position of animals and as a proxy of food-chain length, because consumer values manifest significant trophic-step fractionation in ¹⁵N (2–4‰; 3.4‰ heavier on average than those of their prey) [30]. Finally, a δ¹³C–δ¹⁵N dual-isotope approach allows the understanding of conclusive resource partitioning [26] and provides realistic information on the actual source of the diets that was assimilated into aquatic consumers [31,32,33].

We investigated the spatiotemporal variation of a phytoplankton community in a shallow temperate coastal embayment of Yeoja Bay in South Korea and characterized the community patterns in association with environmental conditions. We further examined the food sources of the filter-feeding ark clam Anadara kagoshimensis (also named A. subcrenata and Scapharca subcrenata), which is extensively cultivated in the shallow (water depth, <10 m) subtidal bottom sediment of the bay. Its massive cultivation in the bay may serve as an important ecosystem engineer that governs lower-trophic-level components, major food-web processes, and material cycling [34,35,36,37]. The main purpose of the present study was to identify shifts in the pattern of resource utilization of A. kagoshimensis in response to changes in the community composition of phytoplankton. Based on distinct isotopic compositions between primary producers (e.g., phytoplankton, microphytobenthos (MPB), macrophyte fragments, and terrestrial/river-borne particulate matter [23,38]), we hypothesized that if the clams use varying degrees of diets in accordance with changes in the phytoplankton community, their tissues will register a consistent shift in isotopic values. To test this hypothesis, we patterned the phytoplankton community across seasonal environmental gradients and examined the clam isotopic values associated with the co-occurring phytoplankton assemblage clusters.

2. Materials and Methods

2.1. Site Description and Sampling Protocol

This study was carried out at two subtidal sites in Yeoja Bay on the southern coast of the Korean Peninsula (Figure 1). The bay is a shallow (mean water depth, 5.4 m) temperate coastal embayment with a semidiurnal tide, with tidal amplitudes of a maximum of 3.5 m at spring tides and a minimum of 1.5 m at neap tides. The bay has a total area of about 320 km² and a north–south length of 30 km and an east–west width of 7.2–21.6 km [39]. The surface sediment of the subtidal zone of the bay is composed of finer mud (silty clay, clayed silt, and clay) with a mean particle size of 8.45 Ø [40]. In the northern part of the bay, the uppermost intertidal zone develops a broad (approximately 10 km²) wetland covered by the common reed (Phragmites australis), and three main streams (the Beolgyo, Isa, and Dong streams) empty freshwater into the bay (approximately 1.6 × 10⁵ m³ d⁻¹). The freshwater discharge is concentrated during the summer monsoon rainy season (June–August) [41]. Moreover, marine water intrudes into the bay from offshore through the southern bay mouth, but water exchange is relatively limited because of the narrow channel at the bay mouth.

This broad subtidal area has been developed for the cultivation of the ark clam Anadara kagoshimensis, which accounts for 90% of the domestic production in Korea [42]. This ark shell culture is based on spat seeds collected in the spawning ground within the bay during summer (July–August) and particulate matter that occurs in nature with no supply of food. Therefore, the spatial and temporal variability of food resources may drive the annual culture yields and the reproductive activity that is responsible for spat production and survival [13].

Based on the spatial and temporal variations in the phytoplankton community and the clam growth dynamics, we chose two sampling stations: Station 1 (S1) has a water depth of 2–5 m and Station 2 (S2) has a water depth of 6–9 m during a low and high tide. Phytoplankton, suspended particulate matter (SPM), sedimentary organic matter (sedimentary OM), and clams were sampled monthly from June 2017 to March 2019, with the exception of December 2019, at the two stations. The samplings were conducted at the end of each month.

2.2. Collection and Treatment of Samples, and Analyses of Environmental Parameters

On each sampling occasion, water temperature and salinity were measured onboard using a multiparameter water quality meter (YSI Pro Plus, YSI Inc., Yellow Springs, OH, USA). Approximately 20 L of water was collected at 1 m below the surface using a van Dorn water sampler, followed by prefiltering through a 180 μm mesh screen to remove zooplankton and any large particles. The collected water was placed in acid-washed plastic bottles, kept on ice in the dark, and immediately transported to the laboratory.

In the laboratory, for measurements of chlorophyll a (Chla) and other photosynthetic pigments, 1 L of water was filtered through a pre-combusted Whatman GF/F filter (GE Healthcare, Buckinghamshire, UK) (47 mm diameter; 0.7 μm pore size) and kept frozen at −80 °C. For SPM measurements, 1 L of water was filtered onto a pre-weighed and pre-combusted Whatman GF/F filter (identical to the procedure used for Chla measurement). The filters were oven-dried at 60 °C for 72 h, and the SPM concentration was calculated based on the difference in weight between before and after drying. The filters containing SPM were combusted at 450 °C for 4 h and weighed again after cooling to room temperature (25 °C) in a vacuum desiccator. Subsequently, the suspended particulate organic matter (SPOM) concentration was calculated based on the difference in weight between before and after the combustion of the filters. On each sampling occasion, duplicate measurements for Chla (and photosynthetic pigments) and SPM at each station were conducted and mean values of individual parameters were presented in the present study.

For measurements of dissolved inorganic nutrient concentration, approximately 40 mL of water was subsampled and filtered in situ through a syringe filter (25 mm diameter; 0.2 μm pore size, cellulose acetate; Advantec, Tokyo, Japan), and the filtrates were kept frozen (−80 °C) until further analyses. Before the nutrient analysis, the water samples were thawed overnight at low temperature (2 °C) and then brought to room temperature (25 °C). Duplicate samples for nutrient measurements were analyzed on a QuAAtro nutrient auto-analyzer (Seal Analytical Ltd., Southampton, UK) according to the procedures developed for phosphate (PO₄³⁻ [43]), ammonium (NH₄⁺ [44]), nitrite (NO₂⁻), nitrate (NO₃⁻), and silicate (Si(OH)₄ [45]).

2.3. Pigment Analysis and Chemotaxonomic Identification of Phytoplankton

Filters for photosynthetic pigment analysis were extracted with 5 mL of 95% methanol for 12 h in the dark at −20 °C and sonicated for 10 min, to disrupt cells. The filters were then ground with a homogenizer and the solution was extracted at 3000× g for 10 min, to remove cellular and glass filter debris. The individual extracts were filtered through a 0.45 mm polytetrafluoroethylene syringe filter, and an aliquot (1 mL) of the filtered solution was mixed with 300 µL of water in an analytical vial. The vial was transferred to an automated sampler. A 100-μL sample of this solution was analyzed using a reverse-phase high-performance liquid chromatography instrument (LC-20AD, Shimadzu, Kyoto, Japan) equipped with a Waters Symmetry C8 column (4.6 × 150 mm; particle size, 3.5 μm; pore size, 100 Å; Waters, Milford, MA, USA) according to the procedure described in [46]. All diagnostic pigments (Chla, chlorophyll b, fucoxanthin, peridinin, 19-hexanoyloxyfucoxanthin, 19-butanoyloxyfucoxanthin, prasinoxanthin, violaxanthin, neoxanthin, diadinoxanthin, alloxanthin, zeaxanthin, lutein, and β-apo-carotenal) were quantified using a spectrophotometer with known specific extinction coefficients [47]. Additional details of the procedures used for the identification of sample peaks and the quantification of individual pigments are provided elsewhere [48].

Based on the initial estimate of the ratios of the diagnostic biomarker pigments, the CHEMTAX program (https://data.aad.gov.au/metadata/records/CHEMTAX, accessed on 1 May 2022) was used to estimate the contribution of eight phytoplankton classes (diatoms, cyanobacteria, chlorophytes, prymnesiophytes, pelagophytes, cryptophytes, dinoflagellates, and prasinophytes) to the total Chla concentration [49,50]. The class-specific input pigment ratios reported previously for various phytoplankton species from Korean waters [51] were used in the CHEMTAX calculation procedure.

2.4. Self-Organizing Map Analysis and Clustering of the Phytoplankton Communities

A Kohonen self-organizing map (SOM) was used to classify the seasonal and spatial patterns of the phytoplankton communities. The SOM performs a nonlinear projection of multivariate datasets onto a two-dimensional space. The SOM analysis employs a competitive and unsupervised artificial neural network algorithm [52,53] and can shape complex community data by training the input data and mapping the clusters of communities on a two-dimensional grid [54,55]. The input datasets on taxa abundances were logarithmically transformed with a base of 10 [log(abundance + 1)] prior to the SOM training, to reduce the variation and skewness of the abundance. The input neurons of the input layer of the SOM are assigned to the output neurons of the output layer consisting of adjacent neurons with neighborhood relations, and the relationship is weighted [53]. Starting from the virtual sites of the input datasets, the random sequences are updated using an iterative training algorithm, to determine the best match with the input neurons, and the further aligning phases are repeated to adjust the neighboring units. Using iterative learning procedures, the best-matching unit and its topological neighbors are finally linked to the corresponding unit of the map. The size of the map grids was confirmed based on two criteria, i.e., (1) Vesanto’s heuristic rule of 5$\sqrt{\mathrm{N}}$ [56], where N is the size of the dataset, and (2) the minimum values of the quantization error (QE) and topographic error (TE) for evaluating the quality of the models [55]. The biomass of each phytoplankton group was visualized based on the gradient distribution of color on the SOM, which represents the role of the species in the establishment of the topology classes. We then applied Ward’s minimum variance method to cluster the resultant SOM units using the Euclidean distance. Finally, a comprehensive understanding of the implied attribution of the environmental variables that defined each community cluster was achieved. The SOM training and cluster analysis were implemented using SOM Toolbox [57] embedded in MATLAB software (Version 6.1, MathWorks, Natick, MA, USA).

2.5. Stable Isotope Analyses

Ark shells and feasible nutritional sources for measurements of carbon and nitrogen stable isotope ratios (expressed as δ¹³C and δ¹⁵N, respectively) were collected on each sampling occasion. For measurements of SPOM stable isotopes, approximately 20 L of water was additionally collected at 1 m below the surface using a van Dorn water sampler (Wildco Instruments, Yulee, FL, USA) and prefiltered through a 180 μm mesh screen to remove zooplankton and any large particles. The water samples were then filtered on a Whatman GF/F filter (25 mm diameter; 0.7 μm pore size). The filtered samples were decarbonated by fuming with HCl overnight in a desiccator, and then dried at 60 °C for 72 h. To compare isotope values with those of SPOM, microphytoplankton (composed mainly of diatoms and dinoflagellates) were collected by obliquely towing a conical plankton net with a 45 cm diameter mouth and a 20 μm mesh aperture at each station. After filtration using a 180 μm mesh screen, the collected microphytoplankton was concentrated in 50-mL test tubes; subsequently, several drops of 0.12 N HCl were added to the samples used for δ¹³C measurements until the bubbling stopped, whereas those used for δ¹⁵N measurements were left untreated. Sediment samples were collected using a 1.2-m² van Veen grab and transported to the laboratory after scraping the top 0.5 cm of the sediment. Sedimentary OM was prepared by adding several drops of 1 N HCl (to remove carbonates) until the bubbling stopped. The sediment samples were then dried in an oven at 60 °C for 72 h and pulverized to a fine powder using a ball mill (Retsch MM200 Mixer Mill, Hyland Scientific, WA, USA).

The ark shells were collected using a clam-harvesting dredge at each station, and their isotopic values were determined for at least three individuals each month. The collected clams were carefully dissected into shells and soft tissue using a stainless knife. Microphytoplankton and ark shell tissues were lyophilized and then ground to a fine powder using a ball mill (Retsch MM200 Mixer Mill, Hyland Scientific, WA, USA).

Filtered samples were packed into tin disks, and the powdered samples were loaded into tin capsules. The δ¹³C and δ¹⁵N values of the sealed samples were measured using a continuous-flow isotope ratio mass spectrometer (Isoprime, GV Instrument, Manchester, UK) interfacing with an elemental analyzer (Eurovector 3000 Series, EuroVector, Milan, Italy). Isotopic values were expressed in delta (δ) notation as parts per thousand (‰) differences from the conventional standard materials (i.e., Pee Dee Belemnite for carbon and atmospheric N₂ for nitrogen), as follows:

$$\delta X=\left[\left(R_{sample}/R_{standard}\right)-1\right]\times {10}^3$$

where $X$ is the ${}^{13}C$ or ${}^{15}N$ value and $R$ is the ratio of ${}^{13}C/{}^{12}C$ or ${}^{15}N/{}^{14}N$. The analytical precision was within 0.15‰ for δ¹³C values and 0.2‰ for δ¹⁵N values based on replicate measurements of reference materials as internal standards calibrated against the International Atomic Energy Agency standards CH-6 (sucrose) and IAEA-N1 (ammonium sulfate).

In the present study, we used the δ¹³C and δ¹⁵N data for the common reed (Phragmites australis) and MPB, which were recently measured in the intertidal flat adjacent to our sampling stations in the bay (unpublished data, C. Kim; c.f. [58]).

2.6. Statistical Analyses

Prior to statistical analyses, the assumption of the normality and homoscedasticity of the variance of all data was checked using the Shapiro–Wilk procedure and Levene’s test, respectively. The significance of differences in environmental variables between sites was tested using a paired t-test or a Wilcoxon test and using a Kruskal–Wallis test among plankton community clusters, followed by a Mann–Whitney U test. One-way analysis of variance (ANOVA) followed by Tukey’s post hoc test was employed to test the significance of differences in the δ¹³C and δ¹⁵N values of organic matter sources and phytoplankton community clusters. These statistical tests were performed using IBM SPSS Statistics (version 21.0, IBM Corp., Armonk, NY, USA).

Differences between stations and among clusters of the phytoplankton community compositions were tested using a permutational multivariate analysis of variance (PERMANOVA, [59,60]) based on the Bray–Curtis similarity. To test for heterogeneous dispersion of the log-transformed data, we performed a permutational analysis of multivariate dispersions. A PERMANOVA pairwise test was used to test differences between each cluster based on p-values using Monte Carlo permutations. All analyses were performed using the PRIMER + PERMANOVA software ver. 6 [59,61].

We performed an indicator species analysis to define the indicator species, which plays a role in the characterization of the classification schemes of given SOM clusters [62,63]. Indicator species are biological indicators of topological groups at sites, which represented clusters in this study. An indicator value (IndVal) index is calculated using the specificity (A) and fidelity (B) of the species between groups and the group combinations [63]. From this calculation, a square root of IndVal is returned instead of the original index. Indicator species analyses were performed using the “indicpecies” package (version 1.7.6 [64]) in R (version 4.1.3 [65]).

2.7. Isotope Mixing Model

The relative contributions of benthic and pelagic sources of organic matter (SPOM, MPB, and Phragmites) to the clam diets in individual clusters of the phytoplankton community were estimated using a hierarchical Bayesian stable isotope mixing model (MixSIAR package, v. 3.1.9) [66]. Since SPOM include both micro- and smaller-phytoplankton and SOM mixes diverse soures of organic matter, we ran a three-source mixing model using pelagic and benthic sources of oganic matter from clusters A, B, and C-D. We ran the mixing model using both isotope values of clams and organic matter sources with cluster as a fixed factor. During the mixing model calculations, we used the isotopic discrimination factors (trophic enrichment factors, TEFs) of 0.8‰ ± 0.2‰ (mean ± SD) for δ¹³C and 3.5‰ ± 0.1‰ for δ¹⁵N between clams and their feasible dietary items, as known for infaunal suspension-feeding bivalves [67]. The model ran Markov chain Monte Carlo sampling based on the following parameters: chain length = 100,000, burn = 50,000, thin = 50, and chains = 3. Diagnostic tests (Gelmin–Rubin, Heidelberger–Welch, and Geweke) and trace plots were applied to examine for any model convergence. The estimated proportions of diet were reported as median (%) and associated 95% credible intervals.

3. Results

3.1. Environmental Conditions

The water temperature in Yeoja Bay was lowest (0.9 and 1.4 °C, respectively) in January 2018 and highest (30.3 and 30.0 °C, respectively) in July 2018 at S1 and S2, representing a seasonal cycle that is typical of the temperate zone (Figure 2A). The salinity range was 28.5–34.7 psu at S1 and 28.4–34.7 psu at S2. There were no significant differences in temperature (paired t-test, t = −0.037, p = 0.970) and salinity (t = −0.753, p = 0.460) between the sites (Figure 2B). The SPM and SPOM concentrations remained constant at ~20 and ~8 mg L⁻¹, respectively, at both stations, with the exception of the abnormal peaks observed in May 2018 and January 2019 (Figure 2C,D). Although a statistically significant difference in SPM concentrations was recorded between stations (Wilcoxon test, Z = −3.042, p < 0.05), the SPOM concentrations showed no significant differences between S1 and S2 (Z = −1.269, p = 0.204). The monthly dissolved inorganic nutrient (NH₄, NO₂, NO₃, PO₄, and SiO₂) concentrations fluctuated irregularly, with no clear seasonality (Figure 2E–I). A Wilcoxon test revealed an absence of significant differences in NH₄, NO₃, PO₄, and SiO₂ concentrations between the stations (Z = −1.154, −0.283, −1.655, and −1.720 for S1; p = 0.248, 0.777, 0.098, and 0.085 for S2, respectively). A significant difference in NO₂ concentration between stations was detected (Wilcoxon test, Z = −2.396, p < 0.05); however, the concentrations were negligible (<0.7 μM) compared with those of NH₄ and NO₃ recorded at both stations.

3.2. Chla Concentration and Phytoplankton Community Composition

The monthly Chla concentration peaked in summer (>5 μg L⁻¹ in June–August) 2017 and fluctuated irregularly in 2018 and 2019 at both stations (Figure 3A,B). The peaks in Chla concentration were followed by a sudden drop to <2 μg L⁻¹. A Wilcoxon test revealed an absence of significant differences in Chla concentration between the stations (Z = −1.756, p = 0.079).

Diatoms predominated at both stations, with an average contribution of 74.0% at S1 and 73.6% at S2. In spring (May–June), the contribution of diatoms dropped to 36.4% at S1, followed by cryptophytes (28.3%) and prasinophytes (18.0%). During that period, the diatom contribution decreased to 26.2% and was replaced by cryptophytes (32.7%) and prasinophytes (20.7%) at S2. The contribution of prasinophytes peaked in September 2018 (23.0%) at S1 and in June 2018 (27.1%) at S2. Sharp increases in cyanobacteria were observed in the spring and fall, with contributions of 12.6% (June 2017) and 14.2% (November 2017) at S1 and 12.9% (June 2018) and 14.2% (October 2018) at S2. Chlorophytes showed high contributions in the fall of 2017 (October, 23.4% at S1; September, 11.3% at S2) and the spring of 2018 (12.8% at S1 and 17.3% at S2). The contribution of dinoflagellates was highest in July (10.0% and 15.0% for S1 and S2, respectively) and December (13.5% at S1 and 12.0% at S2) 2017. Pelagophytes and prymnesiophytes were minor contributors with an average contribution <1% at both stations. Finally, one-way PERMANOVA revealed an absence of significant differences in the phytoplankton community composition between stations (Pseudo-F = 0.359, p = 0.831).

3.3. Patterning the Phytoplankton Communities

After training 42 individual datasets (2 stations × 21 months) of photosynthetic pigment samples, the SOM was generated by a 7 × 4 rectangular grid (28 units, based on a QE of 0.10 and a TE of 0.00) (Figure 4A). The map configuration revealed a clear temporal partitioning of unit samples: the non-bloom season of phytoplankton in the upper left part and the bloom season in the middle and bottom parts of the map. A hierarchical cluster analysis revealed that the SOM units were divided into four clusters A, B, C, and D (Figure 4B).

A PERMANOVA test showed significant differences in phytoplankton community compositions among the clusters (pseudo-F_3,38 = 20.02, p < 0.001), and subsequent pair-wise tests revealed significant differences between the pairs of the clusters (p < 0.05) (

Table 1. Results of PERMANOVA and pairwise tests on phytoplankton community composition of clusters.

Source	df	MS	Pseudo-F	p (perm)
Cluster	3	8537.8	20.02	<0.001
Residual	38	426.51
Total	41
Groups	t		p (perm)
A vs. B	4.271		<0.001
A vs. C	5.427		<0.001
A vs. D	4.072		<0.001
B vs. C	5.152		<0.001
B vs. D	2.523		0.006
C vs. D	2.951		<0.001

). Cluster members indicated that those differences conformed to seasonal successions, not spatial patterns (Figure 4A). All the samples in September–November 2017 belonged to cluster A. Most samples in May–June 2017 and 2018 fell in cluster B. In January–April 2018 and August of both years, most samples were associated with cluster C. In September–November 2018, cluster C was replaced by cluster B with the exception of samples in September at S1 and November at S2. Cluster A reappeared as a dominant cluster in January–March 2019. In July–December 2017, all the samples were covered by cluster D and samples in July 2018 belonged to cluster A.

As shown by the spatial distribution gradients of phytoplankton groups on the SOM units (Figure 4C), a Kruskal–Wallis test revealed significant differences (p < 0.05) among clusters in all the phytoplankton groups (Figure 5). The phytoplankton groups had the lowest concentration in cluster A and the highest in cluster B (Mann–Whitney U test, p < 0.05), with the exception of diatoms and dinoflagellates that had the highest concentrations in clusters C–D and D, respectively. On the SOM, prymnesiophytes, chlorophytes, and prasinophytes exhibited maximum concentrations in the upper-right part of the map, occupying the cluster B units. In turn, cyanobacteria, cryptophytes, and pelagophytes occurred at high values in the mid-to-lower–right part of the map, constituting the main components of cluster B and D units.

IndVal index values revealed the phytoplankton groups that characterized each cluster (

Table 2. Results of the indicator species analysis for the SOM clusters.

Cluster	Group	IndVal	p
Cluster B	Prymnesiophytes	0.86	<0.001
	Chlorophytes	0.82	0.002
	Prasinophytes	0.82	0.002
Cluster C	Diatoms	0.59	0.049
Cluster D	Dinoflagellates	0.90	<0.001

). Prymnesiophytes, chlorophytes, and prasinophytes were significantly associated with cluster B (0.82–0.86, p < 0.01), whereas diatoms were indicative of cluster C (0.59, p < 0.05) and dinoflagellates were an indicator taxa that characterized cluster D (0.90, p < 0.001).

3.4. Characterization of the Environmental Conditions of Clusters

A Kruskal–Wallis test revealed an absence of significant differences (p > 0.05) in temperature, salinity, SPM, SPOM, NO₂, NO₃, and PO₄ among the clusters; by contrast, a significant difference (p < 0.001) in Chla, NH₄, and SiO₂ values was detected among the clusters (Figure 6). A subsequent Mann–Whitney U test revealed that cluster A had the lowest Chla concentration (1.2 μg L⁻¹) compared with other clusters (median: 5.6, 4.5, and 7.2 μg L⁻¹ for clusters B, C, and D, respectively). Moreover, the NH₄ concentration was high in clusters A, B, and C (median: 1.1, 2.6, and 2.0 μM, respectively) compared with that observed in cluster D (median, 0.2 μM; Mann–Whitney U test, p < 0.05). Finally, the SiO₂ concentrations were extremely high in cluster B (median, 29.3 μM) compared with the remaining clusters (median, 3.3–9.1 μM).

3.5. Isotope Signatures of Feasible Diet Components and Clams

ANOVA testing revealed that the δ¹³C and δ¹⁵N values of the primary sources of organic matter differed significantly among them (F_2,89 = 46.216, p < 0.001 and F_2,89 = 53.828, p < 0.001, respectively;

Table 3. δ¹³C and δ¹⁵N values (‰) of primary organic matter sources. The significance of differences in δ¹³C and δ¹⁵N values among potential food sources and SOM clusters was tested using a one-way analysis of variance (ANOVA) followed by the Tukey post hoc test. Data are the mean ± 1SD and the number of samples (N) is in parentheses. PHY, microphytoplankton; SPOM, suspended particulate organic matter; Sedimentary OM, sedimentary organic matter; MPB, microphytobenthos. The same superscript marked with mean ± 1SD values indicates no significant differences in mean values (p > 0.05). * Isotope data for MPB and Phragmites (unpublished, C. Kim).

Variable	F	df	p	Post Hoc Test
				Comparison	p
δ13C	46.216	2, 89	<0.001	PHY vs. SPOM	0.206
				PHY vs. Sedimentary OM	<0.001
				SPOM vs. Sedimentary OM	<0.001
δ15N	53.282	2, 89	<0.001	PHY vs. SPOM	<0.001
				PHY vs. Sedimentary OM	<0.001
				SPOM vs. Sedimentary OM	0.848
	A		B	C & D	Mean ± SD
δ13C
PHY	−21.9 ± 1.4 (12)		−21.3 ± 2.1 (5)	−21.0 ± 1.9 (15)	−21.0 ± 1.7 (32)
SPOM	−21.6 ± 1.3 a (15)		−23.3 ± 0.9 b (5)	−21.1 ± 1.0 a (18)	−21.6 ± 1.3 (38)
Sedimentary OM	−18.4 ± 0.4 (5)		−17.6 ± 0.9 (3)	−18.0 ± 0.8 (12)	−18.0 ± 0.7 (20)
MPB					−13.5 ± 0.7 (12) *
Phragmites					−27.7 ± 0.4 (12) *
δ15N
PHY	9.7 ± 1.3 (12)		10.5 ± 0.8 (5)	9.8 ± 0.7 (15)	9.9 ± 1.0 (12)
SPOM	6.8 ± 1.7 (15)		6.8 ± 1.5 (5)	7.4 ± 1.1 (18)	7.1 ± 1.4 (38)
Sedimentary OM	7.1 ± 1.1 (5)		6.1 ± 0.4 (3)	6.9 ± 1.5 (12)	6.9 ± 1.3 (20)
MPB					9.7 ± 0.6 (12) *
Phragmites					6.2 ± 0.2 (12) *

, Figure 7). Microphytoplankton (−21.0‰ ± 1.7‰) and SPOM (−21.6‰ ± 1.3‰) showed no difference in δ¹³C values (Tukey’s test, p > 0.05), whereas the mean δ¹³C value of sedimentary OM (−18.0‰ ± 0.7‰) was significantly different from those of microphytoplankton and SPOM. The δ¹³C values of these three components (microphytoplankton, SPOM, and sedimentary OM) were intermediate between those of Phragmites (−27.7‰ ± 0.4‰) and MPB (−13.5‰ ± 0.7‰) recorded on the adjacent intertidal flat. The δ¹⁵N values of microphytoplankton (mean, 9.9‰ ± 1.0‰) were clearly different from those of SPOM (7.1‰ ± 1.4‰) and sedimentary OM (6.9‰ ± 1.3‰). The δ¹⁵N values of SPOM, sedimentary OM, Phragmites (6.2‰ ± 0.2‰), and MPB (9.7‰ ± 0.6‰) were consistent among themselves (Tukey’s test, p > 0.05).

An ANOVA further revealed that the δ¹³C values of microphytoplankton and sedimentary OM showed no significant differences among clusters (F_2,29 = 0.108, p = 0.898 and F_2,17 = 1.524, p = 0.246, respectively), whereas the δ¹³C values of SPOM in cluster B (mean, −23.3‰ ± 0.9‰) were slightly lower (F_2,35 = 8.538, p < 0.001) than those of the remaining clusters, i.e., A and C–D (mean, −21.3‰ ± 1.2‰). The mean δ¹⁵N values of the primary sources of organic matter showed no differences among clusters (F_2,29 = 1.337, F_2,17 = 0.547, and F_2,39 = 1.433; p = 0.278, p = 0.438, and p = 0.251 for microphytoplankton, SPOM, and sedimentary OM, respectively).

The ark clam exhibited significant differences in mean δ¹³C values among clusters (ANOVA, F_2,41 = 13.016, p < 0.001) (

Table 4. δ¹³C and δ¹⁵N values (‰) of primary consumer, ark shell A. kagoshimensis. The significance of differences in δ¹³C and δ¹⁵N values among SOM clusters was tested using a one-way analysis of variance (ANOVA) followed by the Tukey post hoc test. Data are the mean ± 1SD and the number of samples (N) is in parentheses. The same superscript marked with mean ± 1SD values indicates no significant differences in mean values (p > 0.05).

Variable	F		df		p
δ13C	13.016		2, 39		< 0.001
δ15N	1.433		2, 39		0.251
	A	B		C & D	Mean ± SD
δ13C
Ark shell	−17.1 ± 0.9 a (17)	−19.1 ±0.4 b (5)		−17.8 ± 0.8 c (20)	−17.6 ± 1.0 (42)
δ15N
Ark shell	11.3 ± 0.6 (17)	10.8 ± 0.3 (5)		11.1 ± 0.7 (20)	11.1 ± 0.6 (42)

). Subsequently, Tukey’s test revealed that the δ¹³C values of clusters were distinct from each other (p < 0.05). The clams of clusters A (mean, −17.1‰ ± 0.9‰) and C–D (−17.8‰ ± 0.8‰) had higher δ¹³C values than did those of cluster B (−19.1‰ ± 0.4‰). Finally, the δ¹⁵N values of ark shells showed no significant differences among clusters (F_2,39 = 1.433, p = 0.251).

3.6. Contributions of Organic Matter Sources to Clam Nutrition

The MixSIAR mixing model calculation revealed that SPOM and MPB were the main contributors to the clam nutrition in all clusters (Figure 8): median, 40% and 52% in cluster A; 29% and 45% in cluster B; and 54% and 39% in clusters C–D, respectively. The contribution of both nutritional sources accounted for a median of 92% and 93% of the clam nutrition in clusters A and C–D, but decreased to 74% in cluster B. This reflected a noticeable descrease in the relative contribution of SPOM to the clam nutrition in cluster B. The contribution of Phragmites was extremely small (8% and 7%, respectively) in clusters A and C–D, but increased substantially up to 27% in cluster B.

4. Discussion

No physical and bigeochemical variables showed spatial variabilities between stations. Water temperature and salinity displayed clear seasonalities but the concentrations of biogeochemical variables (SPM, SPOM, NH₄, NO₂, NO₃, PO₄, and SiO₂) recorded irregular temporal fluctuations. Based on the associations of assemblages, our SOM clustered the monthly phytoplankton community with little spatial variation between stations. Diatoms dominated the phytoplankton community year-round. Based on the SOM and IndVal index, in spring (May–June), the increased abundances of chlorophytes, prasinophytes, and prymnesiophytes replaced the diatom dominance.Such a temporal variation in the phytoplankton community compositions irrespective to the environmental conditions in the bay was surprising. Consistent with the clustering of the phytoplankton community, our isotopic evidence revealed that their contribution to the clam diets varied in association with the community compositions and biomasses of the phytoplankton. The MixSIAR mixing model estimation denoted that SPOM and MPB (consisting mainly of benthic diatoms) served as the main nutritional contributors to the clams. However, a reduction in clam feeding on pelagic and benthic diatoms was observed in accordance with the dominance of small-sized plankton in May–June. Overall results support our initial hypothesis that changes in the community structures of phytoplankton (e.g., dominance of micro- vs. nano- and pico-sized plankton) would shift the availability of preferential phytoplankton taxa (e.g., diatoms [37]) and thereby the consumption of feasible dietary items by clams.

4.1. Phytoplankton Community Patterns

The shallower S1 (~3 m water depth) and deeper S2 (~8 m water depth) stations showed nearly identical fluctuation patterns regarding physical and biogeochemical variables. Water temperature and salinity exhibited the clear seasonal fluctuations that are typical of the Northeast Asian coastal seas, with low salinity in the warm summer monsoon and high salinity in the cold winter season. By contrast, biogeochemical variables were characterized by sporadic peaks of high concentrations (e.g., SPM, NO₃, PO₄, and SiO₂) with no clear seasonality (Figure 2). This result suggests a strong benthic–pelagic link via the effect of external forcing (e.g., tidal currents and wind waves) on the shallow subtidal water, resulting in homogeneity of the water column [68,69]. As indicated by the short residence time (~7 days) of the bay water [41], the water flow through the southern bay-mouth channel may induce a well-mixed water condition in the entire bay system. Nutrient loading via the irregular pulses of freshwater discharge from the dam located in the north of the bay may affect nutrient dynamics in the bay [70].

No remarkable compositional disparities in the community compositions and biomasses of phytoplankton were observed between stations. Rather, our SOM revealed a temporal pattern in the community compositions of phytoplankton. Although the indicator taxa that characterized individual SOM clusters differed between clusters, diatoms were a dominant class in clusters A, C, and D. The fall and winter communities were differentiated into clusters A (September–February) and C (September–November); cluster A was representative of a community during the non-bloom period with the lowest biomass (Chla) in all phytoplankton groups, with no particular indicator taxa, whereas cluster C was characterized by the highest biomass of diatoms (Figure 6). Dinoflagellates served as an indicator class of cluster D, with a transient occurrence, but were detected concurrently with a high abundance of diatoms. Therefore, considering the dominance (and co-occurrence) of large-sized phytoplankton taxa (i.e., diatoms and dinoflagellates, ≥20 μm microplankton) and the prey-size selectivity of clams [24,25], clusters C and D should be merged into a group for further isotopic identification of the clam diets (see later in the Discussion).

Another interesting result was the dramatic compositional shift in phytoplankton communities observed in May–June in both years. Despite the year-round diatom dominance in the phytoplankton community, the dominant occurrence of smaller-sized phytoplankton (0.2–2 μm picoplankton and 2–20 μm nanoplankton), such as chlorophytes, prasinophytes, and prymnesiophytes, in cluster B was surprising and facilitated an early summer bloom. The occurrence of those taxa was associated with warm and highly saline conditions. The extremely high SiO₂ concentrations compared with other clusters likely reflected its reduced consumption caused by the low abundance of diatoms (i.e., siliceous taxa) during that period. The cause of the preference of small-sized phytoplankton for NH₄ nitrogen [71] was unclear in the present study. Because of the absence of conspicuous seasonal and spatial variabilities in the observed nutritional conditions, with the exception of the very few sporadic peaks of inorganic nutrient concentrations mentioned above, it was difficult to define the major environmental determinants of the cluster groups in the present study. One of the most probable explanations for the present clustering of the phytoplankton community may be transport via the renewal of coastal water outside the bay. This interpretation can be evidenced by the indicator taxa of clusters B (small-sized plankton) and D (dinoflagellates), which are the major components of phytoplankton assemblages, especially during spring and summer in the southern coastal waters of Korea [72]. The effect of temporal shifts in the taxonomic (i.e., size) structure of phytoplankton on the dietary changes of clams is discussed in the following section.

4.2. Temporal Variation in the Contribution of Phytoplankton to Clam Nutrition

The δ¹³C values of organic matter components available to the clam nutrition in Yeoja Bay spanned a wide range (−25.9‰ to −15.3‰), in the order of Phragmites, SPOM, phytoplankton, sedimentary OM, and MPB, with their δ¹⁵N values falling within a narrow range (6.9‰–9.9‰). The δ¹³C and δ¹⁵N values of these dietary components measured here were consistent with those reported previously in the southern coastal waters of Korea [23,37,48,58] and were clearly different among components and between seasonal clusters in the case of SPOM. Despite the compositional variation, no significant difference in the δ¹³C values of phytoplankton was observed between community clusters, likely reflecting the persistent dominance of diatoms [23]. The seasonal consistency in the δ¹³C values of sedimentary OM indicates a major contribution of phytoplankton and MPB to the sedimentary OM pool. By contrast, the more negative δ¹³C values of SPOM compared with sedimentary OM, specifically during the period assigned to cluster B, suggest a considerable contribution of ¹³C-depleted organic matter (e.g., small-sized phytoplankton, Phragmites-derived detritus, and terrestrial materials) to the SPOM pool [23,37,58].

The δ¹³C values of clams were arrayed between phytoplankton and MPB, indicating a major contribution of microalgal organic matter to clam nutrition [31]. However, the δ¹³C values of the clams exhibited a shift to more-negative values (mean, −19.1‰) in cluster B of the phytoplankton community compared with clusters A and C–D (mean, −17.1 and −17.8, respectively), indicating a temporal change in their dietary use. Our mixing-model estimation highlighted the fact that phytoplankton and MPB are the main contributors to the clam diets during the period assigned to clusters A and C–D in the phytoplankton community. This result supports the previous finding that benthic diatoms, which are the major constituents of MPB, are also important nutritional contributors to benthic consumers in the bay [31]. In turn, the dietary contributions of SPOM (pelagic microalgal components) decreased, whereas the contribution of Phragmites detritus increased during the period corresponding to cluster B. As mentioned above, our cluster grouping of phytoplankton revealed a prominent temporal shift from micro-sized diatoms in winter to early spring (clusters A and C) to small-sized phytoplankton taxa (i.e., pico- and nano-sized chlorophytes, prasinophytes, and prymnesiophytes) in late spring (May–June, cluster B), when clam spawning and spat recruitment occur. The reduced contribution of SPOM (phytoplankton) to the clam diets during the period of cluster B was linked to the dominance of small-sized taxa.

Both the seasonal community compositions and bloom dynamics of phytoplankton are of particular importance in relation to food availability (both quantity and quality) to clams [73,74]. In Yeoja Bay, the effect of long-term warming on the physiological performance and seasonal cycle of the flesh weight of cultured ark clams has been reported in different manners between winter and summer [13]. Although the improved energy acquisition afforded by the physiological benefits of winter warming promotes the weight gain and gametogenic development of clams, the physiological energy-balance disruption caused by heat stress during summer warming results in progressive weight loss after spawning [13,14]. In combination with the temperature effect, food availability may be crucial in fulfilling clam nutritional requirements [75] during the critical period (late spring to summer), as follows: (1) for gamete development and tissue growth in adults, (2) for post-spawning recovery in adults, and (3) for the survival of spat and juveniles immediately after settlement [13,14], for the persistence of their population size in relation to climate forcing in this warming sea. Our results suggest that, despite the relatively high Chla concentrations, such a dominant occurrence of small-sized phytoplankton during the critical spring period of the clam life cycle may decrease the availability of preferred items (i.e., size-related food quality) and ultimately affect the growth, reproduction, and even survival of the clams.

5. Conclusions

The SOM and subsequent IndVal index analyses of the phytoplankton community revealed a seasonal shift in dominant phytoplankton assemblages in this temperate coastal embayment system. Our isotopic evidence revealed that the trophic importance of pelagic diatoms to clam production may reflect the dominance of these microalgal taxa in the phytoplankton community almost year-round. Benthic diatoms, which are an exclusive taxonomic member of MPB on the broad tidal flat surrounding the bay [76], as well as pelagic diatoms contributed equally to the clam diets. However, when the clams were associated with the dominant occurrence of small-sized phytoplankton taxa (chlorophytes, prasinophytes, and prymnesiophytes) in spring (May–June), the trophic role of diatoms in the clams was diminished and was replaced, to a considerable extent, by other dietary components (i.e., SPOM, sedimentary OM, and Phragmites). This result is consistent with the previous finding that pico- and nano-sized phytoplankton reduce the feeding efficiency of bivalves [19,37]. Our results further indicate that, despite the presence of relatively high Chla concentrations, the dominance of small-sized plankton will lower the food quality for clam feeding. As a result, when linked to a temperature-dependent physiological disruption in clams caused by the warming of the sea [13,14], the decreased food availability might severely threaten clam survival during the critical spring period of the clam life cycle, which possibly explains the recently observed summer mortality [13]. In this context, further research on the combined effects of thermal elevation and lowered food availability on the physiological mechanisms of ark clams [75] is needed for better prediction of the effects of climate change on their growth, reproduction, survival, and population dynamics.