Variability in copepod trophic levels and in feeding selectivity based on stable isotope analysis in Gwangyang Bay off the southern coast of Korea

Trophism (i.e., food resources and trophic levels) of different copepod groups was assessed along a salinity gradient in the temperate estuarine Gwangyang Bay of Korea, based on seasonal investigation of taxonomic results in 2015 and stable isotope analysis incorporating multiple linear regression models. The  13 C and  15 N values of copepods in the bay displayed salinity-associated spatial heterogeneity as well as temperature-related seasonal variations. Both spatial and temporal variations reflected those in isotopic values of food sources. Three major groups (marine calanoids, brackish water 15 calanoids and cyclopoids) had a mean trophic level of 2.2 relative to nanoplankton as the basal food source, similar to the bulk copepod assemblage; however, they had dissimilar food sources based on the different  13 C values. Calanoid isotopic values indicated a mixture of different genera including species with high  15 N values (e.g., Sinocalanus and Labidocera) and relatively low  15 N values (Paracalanus and Acartia). Feeding preferences of different copepods probably explain these seasonal and spatial patterns of the community trophic niche. Bayesian mixing model calculations based on source materials 20 of two size fractions of particulate organic matter (nanoplankton at <20 m vs microplankton at 20–200 m) indicated that Acartia preferred large particles, Paracalanus and Pseudodiaptomus apparently preferred small particles, and Corycaeus was typically omnivorous with low selectivity on particle size. In addition, the carnivorous genus Tortanus predated on copepods without apparent selectivity, Labidocera preferred Acartia to Paracalanus, and Sinocalanus preferred Paracalanus to Acartia and cyclopoids. Overall, our results depict a simple energy flow of the planktonic food web of Gwangyang Bay: 25 from primary producers (nanoplankton) and a mixture of primary producers and herbivores (microplankton), through omnivores (Acartia, Paracalanus, and Corycaeus) and detrivores (Pseudodiaptomus and Euterpina) to carnivores (Tortanus, Labidocera, and Sinocalanus).

field, even though most results of copepod trophism in the MS were similar to previous reports.
Therefore, I would like to recommend the authors to include an additional explanation on a potential limitation which may occur when you apply the current method and assumption to copepod community, in the revised MS.
Response: We thank the reviewer for his/her careful checking of the manuscript and we agree with his/her general assessment of our work and are happy that he found our valuable and positive points in the MS. We agree that the current method has some limitations to study the entire complex planktonic structure. However, the reviewer may misunderstand our data set as we were not just considering the most dominant species only whereas we considered all appearing copepod genera in statistics when we interpreted the dynamics of the entire copepod community or specific taxonomic groups like calanoids and cyclopoids. We believed that copepod genus can be grouped based on similar feeding behaviors and thus the food web structure can be simplified.

Specific Comments:
P6, line 20-22: The author's assumption is questionable. Calanoids consist of many genera or species with various sizes. Even though some large calanoids are not dominant in the sample in terms of abundance, some large calanoids (e.g. Calanus) can P17 line 9-13: The authors used the Bayesian mixing model to estimate the relative contribution of copepods to the carnivore diets, and the prey copepods which were not occurred with predatory copepods according to Table 1 were not considered in the model processing. However, this assumption or process may brings bias when evaluate the prey copepod contribution to predators in reality. The authors did not consider some copepod prey for Labidocera and Sinocalanus, but not Tortanus in Fig. 11. I guess that Labidocera who living on surface also may contact copepods other than Acartia and Paracalanus (for example, according to Table 1, in summer Labidocera rotunda cooccurred with Tortanus as well as Acartia spp.). Therefore, the brackish calanoids and cyclopoid also need to be included in potential prey for Labidocera. The same logic can be applied to Sinocalanus. Although Sinocalanus tellenus dominated in autumn with Paracalanus and Corycaeus, only Acartia was considered as prey for Sinocalanus, but not brackish water calanoid such as Pseudodiaptomus. Please consider all potential prey for Labidocera and Sinocalanus like in the case of Tortanus in Fig. 11A.Also, it is not clear whether the dietary composition of the carnivorous genera in Fig. 10was for a season or for the four seasons. Please specify appropriate season for each carnivorous copepods (e.g. all season or particular season) so that we can guess the potential prey for the carnivorous copepods.
Response: We will try to do so by considering all potential prey for Labidocera and Sinocalanus as suggested by the reviewer. By carefully check the taxonomic dataset, we agree that the brackish calanoids Pseudodiaptomus should be included as a potential food source for Labiocera, as they co-occurred. However, in our taxonomic data set, when we observed Labidocera during the summer, we found that cyclopoid species and didn't occur or may be in extremely low abundance. In such case we don't agree to consider cyclopoid species as potential food source for Labidocera. The dietary composition of the carnivorous genera in Fig. 10 was for all four seasons, as we used the all samples to estimate a mean isotope ratio for each genus. (see Fig. 8A, B, and P17, Lines 11-33, P18, Lines 1-4 and 8-11) P28 Fig.4: Please indicate which genera are the brackish calanoids or marine calanoids in Fig. 5(B) and/or Fig. 5. Also, please specify whether the result of decapods or harpacticoids is for spring and/or winter samples.
Response: I think the reviewer is saying the Fig.5 and Fig.6 (B). We have specified them in figure legend in the revised version. (see Fig. 3  I find a problem in the way the authors estimated the weight differences between cyclopoids and calanoids randomly, as well as assuming that the weight of all calanoid genera was the same. In particular, because the authors have the taxonomic information already, I suggest they do a literature review and obtain the average weight values for each of the copepod genera/species used in the study, and apply these to the bulk regressions. I believe this is especially important as the authors are trying to extrapolate significantly more results than what they measured (i.e. genera-specific isotope values from a mixed community), that the approach be as precise as possible.
Response: We thank the valuable recommendation of the reviewer. We have tried to do a literature review and obtain the average weight values for the important genera used in our study. (see Supplementary Table S1 in the revised ms) In general I appreciate the effort to expand upon simple d13C and d15N bulk measurements for more detailed information on a community. However, in the case of copepods, if the authors do/did intend to investigate these relationships in detail, why not simply measure the values of individual genera? They state that too much material is required, but methodological advances these days mean that an individual Calanus female can indeed be analyzed (60ugC, 10ugN), as 5 ugN is typically the lower limit of standard bulk analyses (and low-N methods methods have been developed to go down to 1 ugN). Cyclopoids would require greater number, but following the authors assumptions of 0.1<x<1, that would be about 10 individuals. When certain problems arise, such as Paracalanus and Sinocalanus having lower d13C values than any measured prey, it would seem the authors acknowledge them, but then continue their analyses, e.g. calculate a TL (presumably based on prey that has been shown to not be consistent with their isotope values) in the same was as for the other genera.
Response: We agree to the reviewer's consideration. However, we missed to do so when the investigation was conducted. As state in the Introduction part, it is too timeconsuming to obtain enough weight for specific genus and isotope analysis for different subgroups requires great expertise in isolating species from highly complex mixtures. Besides, we have tried to simplify the sampling way so that some monitoring departments may follow. (see P6, Lines 13-25 in the revised ms) One gets the sense that by plugging it into GAMs and regression models, the error sources and magnitudes are lost. I would like to see a quantitative test of the biases inherent in this Bayesian model, and how confident the authors can be that this approach is recovering the actual copepod diets. Given this approach and the number of assumptions that lie within, uncertainty relating to the model (as well as replication, independently) should be presented, discussed, and assessed explicitly with the other sources of uncertainty. This should be done with both the particle feeders and the carnivorous species, and the effects of including or excluding different species types should also be assessed.
Response: The reliability of Bayesian mixing model is fully discussed in literature (Phillips and Koch, 2002;Phillips and Gregg, 2003;Moore and Semmens, 2008;Ward et al., 2010;Parnell et al., 2010Parnell et al., , 2013. We are not good at extrapolating this model. However, we have tried to present more about the details we used, e.g. the replication, trophic enrich factor, sources just like the reviewer suggest. (see P7, Lines 11-27 in the revised ms) Finally, consistent with the point I discuss above, the authors mention a 'simple energy flow' in the abstract and discussion. But I wonder if this methodological approach allows for more complex flows. The actual isotopic values were not measured, but inferred from mass balance of dominant genera, and Bayesian approaches, and the violation of the underlying assumptions was not determined. How would a more complex picture emerge? In fact, the problem of Paracalanus and Sinocalanus having lower d13C values could hint at more complexity, yet it is assumed perhaps that this is due to unmeasured food sources and then ignored. I think if the authors address the issues posed above (and specifics below), the MS is suitable for publication.
Response: We admit that a complex picture cannot be fully understood from this paper.
The estimated values only can provide a mean value and a standard error, thus they cannot explicit exactly the same situation of different seasons and no dynamics picture can be found. But the mechanism to estimate the isotopic ratio of a mixing sample, which is mixed by different species with different masses, is clear. Our results only suggest the potential trophic position of those examined genera. For Paracalanus and Sinocalanus, their low 13C was significantly estimated from the samples containing certain amount of them. We believe the data are true and correct. Although we fail to provide information of all potential food items for them, they show the role in interacting with plankton and other copepods. Frankly speaking, we didn't intend to give detail information of the biology of each genus, but aims to investigate their potential roles in regulating the abundance of the two size fractions of plankton. (see P18, Lines 20-424 in the revised ms) Response: Agree to change. (see P1, Lines 14-15 in the revised ms) Introduction P2-0. "With broad feeding spectra and flexible feeding strategies, the bulk copepod assemblage is omnivorous depending on dominant species or group". Omnivorous or what? Consider changing to something like 'displays varying degrees of herbivory/ omnivory/carnivory, depending on dominant: : :' Response: Agree to change. (see P2, Lines 4-5 in the revised ms) P2-5. "In turn, TLs of a diverse: : :". I assume the authors here refer to the average trophic position of the assemblage, and thus should be 'TL' (singular). "Because copepods play a fundamental role in feeding on phytoplankton as primary consumers".
Consider re-phrasing as 'Because copepods rely significantly on phytoplankton as prey', otherwise the expectation of this phrase is that the second half will refer to the top down effect of copepods on phytoplankton, and not the bottom-up effect of phytoplankton on copepods. 'feeding on phytoplankton as primary consumers, so the seasonal and spatial'.
Response: Agree to change. (see P2, Line 9 in the revised ms) P2-15. "Therefore, the assessment of the trophic position (: : :) of copepods within a complex planktonic food web is critical in predicting the ecological relationships between predator and prey". This phrase seems redundant, isn't the study about assessing these ecological relationships? I don't understand the prediction part.
Response: To avoid confuse, we have revised in predicting to to understand. (see P2, Line 21 in the revised ms) P3. 0. "In contrast, the d15N values of primary producers increase from being nutrient sufficient (high fractionation) to nutrient-limiting (low fractionation) and are especially high in anthropogenic wastewater nitrogen inputs". Would the later simple swamp the fractionation effect? The literature on ïA˛d'15N of different nutrients in the ocean (nitrate, ammonia, urea) shows ranges that are much larger than fractionation factors, e.g. these vary by about 20‰ compared to 3.4‰ of fractionation. Can you comment on how much you expect the source to vary along the river gradient?
Response: The parentheses here indicate the consequences. For example, rich nutrients will cause high fractionation of primary producer. And they would continuously accumulate their δ 15 N in the cells. We don't have the data on the variations of δ 15 N of source, as we expect theδ15N of source will vary from 0 to 13‰ based on the variances of POM. (see P3, Line 5-7) Materials and methods P4-25. Could you mention the average volume filtered per tow, as the net was equipped with a flow meter?
Response: Sure, we can provide this. (see P5, Lines 3-4 in the revised ms) P5-5. "water samples were transported to the laboratory as soon as possible". Please give a time estimate.
Response: 12 hours driving. (see P5, Line 14 in the revised ms) P6-5. The analytical precision of 0.2‰ and 0.3‰ for d13C and d15N, respectively, seems a bit high. Could you estimate what is the lowest change in TL that you can estimate based on this instrument error?
Response: We have re-checked the precision of the instrument during the period we measured the samples. The analytical precision should be 0.1 and 0.05 ‰ for δ 13 C and δ 15 N, respectively. Based on the equation of trophic level, the lowest change will be less than 0.02 TL. (see P6, Line 11 in the revised ms) P6.15. The weight difference between cyclopoids and calanoids was generated randomly. I don't understand why the information from the species identification was not used for this purpose. What is the error associated with this type of computation? I would really suggest the authors do a literature search of the mean weights of the difference species and genera enumerated in their samples, and use this information to estimate both cyclopoid/calanoid weights, and the weights of the different calanoid genera. If the composition has already been estimated, it makes no sense to make these assumptions that only introduce greater error into an already indirect way of estimating species stable isotope composition.
Response: As mentioned above, we have doe a literature search. (see Supplementary   Table S1 and P6, Lines 13-22 in the revised ms) P7-20. 'fractionation factors used in the model estimation were calculated from TLs'. I don't understand this statement, it sounds like 3‰ and 0.5 ‰ were assumed (logically) and not calculated. Please clarify.
Response: The sentence of "fractionation factors used in the model estimation were calculated from TLs' means that the fractionation factor between two trophic levels, based on the difference of each δ 15 N value. Or it can also be understood by this way, it is multiplying the difference of two trophic levels (calculated from different δ 15 N value) with 3‰ and 0.5 ‰ for δ 15 N and δ 13 C, respectively. (see P7, Lines 22-23 in the revised ms) P9-15. "Copepod d15N : : : being much more consistent with the pattern of microplankton than that of nanoplankton". This seems true for the summer ïA˛d'15N values, and quite the opposite for the winter values. Regardless, there is such high variability that it is hard to tease out any clear pattern of spatial/seasonal co-variability.
Response: Yes, it is true. We have removed this unsuitable sentence. (see P9,  in the revised ms) P9-20. The GAM result is very interesting. Perhaps it reflects the food-web processes that affect d15N disproportionally and were not included in the GAM?
Response: Here the deviances explained by GAM suggest those factors combined to influence the dependent factor. For δ 15 N, the deviance explained was relatively lower suggesting that other factors which were not included would contribute another 23% of the deviance of δ 15 N. But we don't know what are them. The understanding of the reviewer is right. (see Table 2, Table 3, and P9, Lines 25-P10, Line 32 in the revised ms) P10-20. It is not clear to me how the trophic levels of brackish copepods can be calculated, when their 13C values do not support the sampled nanoplankton and/or microplankton as their food source. I also don't understand how later in figure 6 they show up enriched, but in figure 4 they are depleted with respect to this food source.
The differences between these two figures should be stated clearly as they show different results.
Response: The brackish copepods in this study were defined by empirical taxonomy, including Pseudodiaptmus and Sinocalanus. The trophic levels were calculated by the formula shown in M&M and figure legend. We admit that there may be some confuse, while we don't think the two figures show different results. Firstly, the results of brackish copepods were averaged from Pseudodiaptmus and Sinocalanus, while our result showed more insight information that they were different.
Secondly, as mentioned above, we didn't intend to give detail information of the biology of each genus, but aims to investigate their potential roles in regulating the abundance of the two size fractions of plankton. Thus, we were unable to provide detail information of all potential diet sources of them, whereas they still have enrichment factor on lower trophic levels such as nanoplankton and microplanton.
That's not right, it should say something like 'enrichment values for marine and brackish water: : : feeding on nanoplankton'.
Response: Sorry for a mistake here. We have revised as suggested. (see P11, Lines 6-12 in the revised ms) P11-5. I disagree with the statement (based on the figure) that "the proportions of the two size fractions of POM averaged from all four seasons contributing to copepod diets at different stations were also distinctly different except for station 8 (Fig. 8)". It seems that the error bars overlap at station 1 (hence not different), and stations 6 and 7. I might be missing something but then it should be clarified.
Response: Yes, we found that error bars were indeed overlapping. We have revised this conclusion. (see P11, Lines 19-30 in the revised ms) P11-10. Does 'spring data available' mean 'only spring shown'?
Response: Yes, we obtained enough amounts of decapods for isotopic analysis only at the spring. However, as suggested by another reviewer, we decide to delete the part of decapods as it was not related to the topic of this study. (see P11, Line 27 in the revised ms) The authors discuss size-selective feeding of calanoids in the context of 'filtering efficiency', yet they are not true filter feeders, they are suspension feeders that trap and handle particles (Paffenhofer et al, 1982, Mar Bio 67:2), which has different implications for particle handling. This is an important distinction that should be observed throughout the MS.
Response: OK, we have carefully checked the whole MS and change to feeding efficiency (see P17, Line 7 in the revised ms) Discussion P13-0. It seems to me that the sewage explanation deserves a bit more attention. If the authors can't rule it out it means that this could contribute substantially and swamp the other subtle processes discussed in the 15N-enriched ammonia section.
Response: Yes, we also believe that sewage was important for 15 N accumulating.
However, we didn't have direct data to support our speculation. Thus we have changed the sentence to The input of sewage-derived 15  P13-5. "Furthermore, the fractionation effect of phytoplankton will be reduced when phytoplankton became nitrogen-limited and take up nitrogen with little fractionation". I am unsure that this effect could be significant in a coastal areas such as this one.
Moreover, if phytoplankton reduce their fractionation, it would mean that their 15 N will tend to be higher (as they choose the lighter 14 N), and thus doesn't explain this decreasing trend.
Response: Yes, we agree that nutrient-limiting is not frequently happened in coastal area.
However, substantial reduction of nutrients from different seasons or from different stations and the mis-match of high phytoplankton and low nutrients were normal.
When phytoplankton reduce fractionation, they will select more ligther 14 N in cells thus they will show a reducing ratio of 15 N in cells (Cifuentes et al., 1988;Fogel and Cifuentes, 1993;Granger et al., 2004). To remove such confuse, we have revised this sentence to Furthermore, the fractionation effect of phytoplankton will be reduced when nutrients substantially decreased and phytoplankton would take up nitrogen with little fractionation and stored relatively light of nitrogen isotope. (see P13, Lines 24-25 in the revised ms) P13-10. I would like to see table with the GAM results. It would be nice to have these presented first in the results, and later discussed. It would also be interesting to see the different variables tested and the ones found to be significant within this table.
Response: Agree to revise. We have tried to show the table in Results. And as suggested by another reviewer, we will move the GAM figures ( I am somewhat confused about the discussion of trophic levels of the copepods Paracalanus and Sinocalanus. The authors state that their ïA˛d'13C values are lower than all measured food sources, which would imply that their food source has not been adequately measured. How then are these organisms included in the trophic level (TL) component of the paper? A bit of clarification on this topic would really help the reader.
Response: Agree to do so. The trophic level in this study was defined as trophic position relative to nanoplankton, which was considered as the trophic baseline. (see P10,  in the revised ms) P17-10. This paragraph explaining the Bayesian mixing model methods/results should be moved to the results section.
Response: Agree to revise. Instead, we have added some more discussion on this part.

General Comments
This manuscript provides results from seasonal and spatial variation in the stable isotopes 13C and 15N of POM and copepods along a salinity gradient in Gwangyang Bay, off the southern coast of Korea. The authors combined this information with linear mixing models, Bayesian isotopic mixing models and generalized additive models to derive a statement on food selectivity and trophic level of copepods. In general, this manuscript is very well structured and provides valuable information on the flow of matter through the food web. Still, some concerns have to be clarified before publication.
Response: We appreciated the positive comments of the reviewer and have followed the suggestion to improve the manuscript. Response: Here the "highly mixed species" means the assemblage contained too many different species and those species had similar size. So such species were hard to be sorted out from the assemblage based on current microscopic technique. To remove confuse, we have revised it to high diversity of the assemblage and ….

Specific comments.
(see P3, Line 21 in the revised ms) 3. Page 3, line 21: Instrument sensitivity has increased and compound specific analysis (CSI) of stable isotopes in amino acids make it possible to track diets of mesozooplankton and determine their trophic position.
Response: Yes, of course. We admit that highly developed instrument can do so. But for doing so, researchers still need taxonomic expertise to sort out the species from a complex mixture to prepare the sub-sample. It requires a lot of lab processing Response: Based on revised estimation, we have tried to provide such simplified figures.
(see Fig. 9 and P8, Line 20-24in the revised ms) I hope that these revisions are satisfactory and that the revised version will be acceptable for publication in Biogeosciences.
Sincerely yours, Abstract. Trophic preference (i.e., food resources and trophic levels) of different copepod groups was assessed along a salinity gradient in the temperate estuarine Gwangyang Bay of Korea, based on seasonal investigation of taxonomic results in 2015 and stable isotope analysis incorporating multiple linear regression models. The  13 C and  15 N values of copepods in the bay displayed significant spatial heterogeneity as well as seasonal variations, which were indicated by their significant relationships with salinity and temperature, respectively. Both spatial and temporal variations reflected those in isotopic values of food 15 sources. The major calanoid groups (marine calanoids and brackish water calanoids) had a mean trophic level of 2.2 relative to nanoplankton as the basal food source, similar to the bulk copepod assemblage; however, they had dissimilar food sources based on the different  13 C values. Calanoid isotopic values indicated a mixture of different genera including species with high  15 N values (e.g., Labidocera, Sinocalanus, and Tortanus), moderate values (Calanus sinicus, Centropages, Paracalanus, and Acartia), and relatively low  15 N values (Eurytemora pacifica and Pseudodiaptomus). Feeding preferences of different 20 copepods probably explain these seasonal and spatial patterns of the community trophic niche. Bayesian mixing model calculations based on source materials of two size fractions of particulate organic matter (nanoplankton at < 20 μm vs microplankton at 20-200 μm) indicated that Acartia and Centropages preferred large particles, Paracalanus, Calanus, Eurytemora, and Pseudodiaptomus apparently preferred small particles. Tortanus was typically carnivorous with low selectivity on different copepods. Labidocera preferred marine calanoids Acartia, Centropages, and harpacticoids, and on the 25 other hand, Sinocalanus and Corycaeus preferred brackish calanoids Paracalanus and Pseudodiaptomus. Overall, our results depict a simple energy flow of the planktonic food web of Gwangyang Bay: from primary producers (nanoplankton) and a mixture of primary producers and herbivores (microplankton), through omnivores (Acartia, Calanus, Centropages, and Paracalanus) and detrivores (Pseudodiaptomus, Eurytemora, and harpacticoids) to carnivores (Corycaeus, Tortanus, Labidocera, and Sinocalanus). 30

Introduction
Mesozooplankton constitute essential trophic mediators of marine food webs in transferring energy and materials by linking the microbial food web to higher trophic levels. Copepods are a diverse assemblage dominating mesozooplankton communities.
With broad feeding spectra and flexible feeding strategies, the bulk copepod assemblage displays varying degree of herbivory/omnivory/carnivory, depending on dominant species or groups (Graeve et al., 1994;Sell et al., 2001;Turner, 2004;5 Vadstein et al., 2004;Gifford et al., 2007;Chen et al., 2017). The role of copepods in planktonic food webs can be determined by their overall trophic levels (TLs) relative to primary producers. In turn, the TL of a diverse copepod assemblage is balanced from different groups with different feeding preferences and are ultimately determined by species composition. Because copepods rely significantly on phytoplankton as prey, the seasonal and spatial changes in the composition and availability of phytoplankton determine the abundance and feeding behavior of the copepod assemblages. 10 The most dominant copepod species, such as Neocalanus, Calanus, Temora, and Paracalanus, are filter-feeders that perform a size-selective feeding behavior depending on particles effectively retained by feeding appendages of copepods. Large phytoplankton (> 20 m; mainly diatoms and dinoflagellates) are generally grazed at high rates by copepods, as shown by many field studies in coastal and estuarine waters (e.g., Liu et al., 2005a, b;Chen et al., 2017). Many other field studies have reported that omnivorous species dominate copepod assemblages because of high feeding selectivity on larger 15 microzooplankton that are considered to have higher nutritional quality (e.g., Berk et al., 1977;Fessenden and Cowles, 1994;Calbet and Saiz, 2005;Gifford et al., 2007;Chen et al., 2013). These omnivorous copepods might induce increases in phytoplankton levels indirectly through trophic cascades as they graze intensely on microzooplankton (e.g., ciliates and heterotrophic dinoflagellates) (Nejstgaard et al., 2001;Stibor et al., 2004;Sommer and Sommer, 2006;Zöllner et al., 2009;Chen et al., 2011Chen et al., , 2013. Therefore, the assessment of the trophic position (herbivores, omnivores, or carnivores) of copepods 20 within a complex planktonic food web is critical to understand the ecological relationships between predators and prey.
Stable isotope analysis (SIA) is a reliable technique providing insight into the trophic positions of copepods relative to basal food sources (Grey et al., 2001;Sommer et al., 2005;Hannides et al., 2009;Kürten et al., 2011). Isotopic comparisons with food sources enable us to analyze prey selectivity during predators' feeding history as well as within food web structures (Fry, 2006;Layman et al., 2012). In general, the carbon stable isotope ratio ( 13 C) can be useful for tracing food sources because of 25 small fractionation (0.5-1‰ per TL) during trophic transfer, particularly when different food sources at a given period in a specific system have distinct  13 C values. By contrast, the nitrogen stable isotope ratio ( 15 N) can be useful for estimating relative TLs because  15 N values of consumers generally increase with TL (an average 3.2‰ of enrichment per TL; Post, 2002;Michener and Kaufman, 2007). The development of linear mixing models and the Bayesian mixing model has allowed researchers to predict the proportions of different food sources in the diets assimilated by grazers (Phillips and Koch, 2002;30 Phillips and Gregg, 2003;Moore and Semmens, 2008;Ward et al., 2010;Parnell et al., 2010Parnell et al., , 2013. Coastal and estuarine environments often experience rapid fluctuations of inorganic carbon and nitrogen inputs in response to diverse oceanographic processes (e.g., coastal currents, upwelling, tidal mixing, and river discharges), which drive spatial and seasonal heterogeneities in biogeochemical dynamics and isotopic signatures (Rolff, 2000). Indeed, the  13 C values of suspended particulate organic matter (POM) in estuarine systems increase progressively from the head to the mouth of each estuary because of the lower  13 C values in terrestrial carbon or sewage materials through river discharge (Cifuentes et al., 1988). In contrast, the  15 N values of primary producers increase from being nutrient-sufficient (high fractionation) to nutrientlimiting (low fractionation) and are especially high in anthropogenic wastewater nitrogen inputs (McClelland et al., 1997). In 5 addition, different phytoplankton groups utilize different nitrogen sources with different enrichment factors, possibly offering different isotopic pools to grazers (Gearing et al., 1984;Rolff, 2000;Montoya et al., 2002). For example, diatoms primarily utilize nitrate with varying fractionation factor on 15 N (0.7-6.2‰) depending on species (Waser, et al, 1998;Needoba et al., 2003), while flagellates primarily utilize ammonia with enrichment factor of 6.5-8‰ (Montoya et al., 1991).
Given that isotopic values of copepods vary in association with copepods' food source by one or two increases in TL values, 10 seasonal and spatial patterns generally follow the trends of their food sources or dominant prey (Grey et al., 2001;Montoya et al., 2002;Kürten et al., 2011). Higher  15 N values of copepods caused by fractionation rather than food source or by averaging from mixed food sources are evident considering the lowered isotopic values of fecal pellets (Checkley and Entzeroth, 1985;Checkley and Miller, 1989;Tamelander et al., 2006). Furthermore, the effect of the microbial food web on the elevated  15 N values of copepods cannot be ignored (Rolff, 2000;Kürten et al., 2011). Therefore, variations in isotope signatures of both 15 copepods and POM (including phytoplankton, bacteria, ciliates, and detritus) help to depict the biogeochemical cycles of specific systems (Grey et al., 2001;Montoya et al., 2002;Francis et al., 2011). Nevertheless, because copepods graze preferentially on larger phytoplankton (diatoms and dinoflagellates) and microzooplankton (ciliates and heterotrophic dinoflagellates), we hypothesize that isotopic values of the copepod assemblage will be much closer to those of larger rather than smaller food source plankton. 20 However, high diversity of the assemblage and size overlap among different species make it hard to determine the relative trophic positions of different subgroups or species. Isotope analysis for different subgroups requires great expertise in isolating species from highly complex mixtures. Moreover, the number of individuals of a specific genus is often insufficient for analysis because of limited instrument sensitivity. Thus, to our knowledge, direct comparisons of different mesozooplankton groups or copepod species are seldom found in the literature (Schmidt et al., 2003;Sommer et al., 2005;Hannides et al., 2009). Here, 25 we estimated isotope values of different copepods by mass balancing linear mixing models from values of bulk samples and taxonomic data of copepods. The allocated masses of calanoids and cyclopoids were achieved from literature and empirical formulas.
Overall, we aimed (1) to understand the seasonal variations and spatial heterogeneity of copepod  13 C and  15 N values in a temperate estuarine system (Gwangyang Bay, Korea); (2) to compare the trophic positions of different copepods; and finally 30 (3) to elucidate the compositions of two major size classes (<20 m and >20 m) of POM in grazer diets. The dietary composition (nano-vs microplankton) of copepods was estimated using Bayesian isotopic mixing models (Parnell et al., 2010(Parnell et al., , 2013. The results of this study will provide insights into trophic preference information (i.e., food resources and trophic levels) for different copepod groups and help in understanding the biogeochemistry of this estuarine system.

Study area
Gwangyang Bay is a semi-enclosed bay system, located on the southern coast of the Korean Peninsula, and is one of the most 5 industrialized coastal areas exposed to anthropogenic pressure. It starts from Seomjin River through Yeosu Channel (between Yeosu Peninsula and Namhae Island) to open ocean (the East China Sea). The bay area covers approximately 145 km 2 and water depth is generally shallow at 2.4-8.0 m in the northern upper-middle Seomjin River estuarine channel compared with 10-30 m in the deep bay channel (Kim et al., 2014). The annual freshwater discharges of Seomjin River are 10.7-39.3  10 8 tons. The seasonality of nutrient input from the catchment area (ca. 5  10 3 km 2 ), including agricultural and forested land, is 10 profound (Kwon et al., 2002). The wet season starts from late spring and the discharge peaks during the summer monsoon period.
Accordingly, the maximum median river discharge varies from 30-95 m 3 s −1 in the dry season to 300-400 m 3 s −1 in the summer monsoon, with an annual mean of c.120 m 3 s −1 (Kim et al., 2014). The tidal cycle of the bay is semidiurnal with maximum ranges of 3.40 m during spring tides and 1.10 m during neap tides. Tidal currents from the Yeosu Channel also strongly 15 influence the system and approximately 82% of the Seomjin River flux is discharged toward this channel. Overall, increasing industrial pollution facilitates eutrophic conditions in the estuarine and related bay waters. Diatoms dominate in the phytoplankton community and density is high in the middle part of the Gwangyang Bay (Kim et al., 2009; data from our parallel study not shown). The distribution patterns of copepods in the Seomjin River during summer were represented by three main salinity zones: an oligohaline zone (predominated by Pseudodiaptomas koreanus, Sinocalanus tenellus, and 20 Tortanus dextrilobatus), a mesohaline zone (predominated by Acartia ohtsukai and Acartia forticrusa), and a polyhaline zone (predominated by Acartia erythraea, Calanus sinicus, Centropages dorsispinatus, Labidocera rotunda, and Paracalanus parvus) (Park et al., 2015).

Sampling and processing
Surface water and net-tow samples were collected seasonally (February, May, August, and November) at nine stations from 25 the head to the mouth of Gwangyang Bay in 2015 (Fig. 1). The stations were chosen based on salinity regime and different geographic characteristics. Stations 1-3 were located in the Seomjin River, stations 4 and 5 were in Gwangyang Bay (the middle part of the estuary), and stations 6-9 were located from the offshore deep-bay channel to the southern mouth of the estuary. On each sampling occasion, water temperature and salinity were determined in situ using an YSI Model 85 probe (YSI Inc., Yellow Springs, OH, USA). 30 Zooplankton taxonomic samples were collected during daytime by net towing using a plankton net (45 cm diameter, 200 m mesh size) equipped with a flowmeter (Model 2030R Mechanical Flowmeter, General Oceanics Inc., Miami, FL) and gently hauled horizontally at a subsurface depth of 0.5-1 m with the ship speed at about 1 knot (0.5 m s -1 ). The average volume filtered per tow was 16.7 ± 5.1 m 3 (mean ± se). Samples were fixed in formalin solution with a final concentration of 5% and then identified and enumerated under a stereomicroscope (SMZ 645; Nikon, Tokyo, Japan) in the laboratory. At each station, 5 one additional net tow was collected for isotope analysis. After collection, specimens were transferred immediately into plastic bottles and preserved in a refrigerator (4 °C) until analysis. In the laboratory, subsamples were picked out from the mixed zooplankton samples under a dissecting microscope. Easily distinguishable zooplankton groups such as harpacticoids were separated from a mixture of calanoids and cyclopoids. Copepodites were counted but grouped together with adults. All subsamples were lyophilized and then homogenized by pulverizing them with a mortar and pestle before isotope analysis. 10 POM in surface water (0.5-1 m depth) was collected using a 5-l Niskin bottle at a midday high tide at the same time as zooplankton collection. Approximately 20 l seawater collected was first screened through a 200 m Nitex mesh to remove zooplankton and large-sized particles. The pre-screened water samples were transported to the laboratory as soon as possible within 1-2 h. In the laboratory, water samples were filtered again through a 20 m Nitex mesh and then filtered onto precombusted (450 °C for 4 h) Whatman GF/F glass fiber filters to determine isotope ratios of fine POM (< 20 m) representing 15 pico-and nano-sized plankton. To obtain enough plankton cells for isotope analysis of coarse POM ( 20m), we collected POM samples by net towing with a plankton net of 50 cm diameter and 20 m mesh size. After collection, each sample was pre-filtered through a 200 m Nitex mesh to remove large particles and zooplankton. Both size fractions of samples were prepared in duplicate. Samples for  13 C measurements were acidified by fuming for about 5 h over concentrated HCl in a vacuum desiccator to remove carbonates, while the samples for  15 N measurements were not acidified. All the samples were 20 lyophilized and pulverized with a mortar and pestle before isotope analysis.
For chlorophyll a (Chl a) determination, 1-l subsamples of surface water were filtered through Waterman GF/F glass fiber filters. The filters for Chl a (including other photosynthetic pigments) were extracted with 95% methanol (5 ml) for 12 h in the dark at −20 °C and sonicated for 5 min to foster cell disruption. Aliquots of 1 ml of the supernatants were mixed with 300 l of water; 100 l of this solution was analyzed by reverse-phase high performance liquid chromatography (HPLC, LC-20A 25 HPLC system, Shimadzu Co., Kyoto, Japan) using a Water Symmetry C8 (4.6  150 mm, particle size: 3.5 m, 100 Å pore size) column (Waters, Milford, MA, USA) and a method derived from Zapata et al. (2000). Quantification of standard pigments was calculated by spectrophotometer with the known specific extinction coefficients after Jeffrey at al. (1997). Sample peaks were identified based on their retention time compared with those of pure standards. Further details on analysis, calibration, and quantification have been given elsewhere (Lee et al., 2011;Kwak et al., 2017). 30

Isotope analysis
For measurements of carbon and nitrogen stable isotope ratios, all pre-treated samples were analyzed using a continuous-flow isotope ratio mass spectrometer (CF-IRMS; Isoprime100, Cheadle, UK) connected to an elemental analyzer (vario Micro cube, Hanau, Germany) following the procedure described by Park et al. (2016). Briefly, powdered samples were sealed in tin combustion cups and filter samples were wrapped with a tin plate. All prepared samples were put into the elemental analyzer to oxidize at high temperature (1030 °C). CO2 and N2gaseswere introduced into the CF-IRMS with the carrier being helium gas. Data of isotope values are shown in terms of X, indicating the relative differences between isotope ratios of the sample 5 and conventional standard reference materials (Vienna Pee Dee Belemnite for carbon, and atmospheric N2 for nitrogen), which were calculated by the following equations:X = [( / ) − 1] × 10 3 , where X is 13 C or 15  analytical precision for 20 replicates of urea were approximately  0.1‰ and  0.05‰ for  13 C and  15 N, respectively.

Data analysis
We used multiple linear regression to assign the isotopic value to each species from a mixture sample. As we did not measure the dry weights each species directly, we calculated them using published empirical equation of the relationship between body length and dry weight of each species and we searched literature for the range of body length and their living environment of 15 each species in the World of Copepods Database and related references ( Assuming that the proportion we picked to do isotope analysis was the same as the proportion in samples for composition analysis, once the weights of different groups or different genera/species were assigned, we computed the isotope ratio for each group by multiple linear regression models as in the following equations: 25 where is the weight of the total community and 1are the weight of different groups or genera of each group.
 13 1 - 13 and  15 1 - 15 are the  13 C and  15 N values of each group or genus, respectively. We used R software to do the estimation using the whole sampling data set. Insignificant results for sparse species, such as Bestiolina coreana, 30 Clausocalanus furcatus, Paraeuchaeta plana, Oithona davisae, and Oncaea venella, found in Gwangyang Bay are not tested separately, while they were incorporated to respective groups based on their living environment.
Given that the isotopic values of consumers come from their diets and thereby from mixed proportions of different sources, the proportions of each source could be simulated by linear mixing models with a fractionation factor (also called a trophic enrichment factor). For instance, a mass balance mixing model is given by: 5  13 = f 1  13 1 + f 2  13 2 + ⋯ + f δ 13 + , f 1 + f 2 + ⋯ + f = 1, where f 1 −f are the proportion of different sources, and and are trophic enrichment factors for  13 C and  15 N values, respectively. 10 Here, a Bayesian isotopic mixing model (available as an open source Stable Isotope Analysis package in R: SIAR) was performed to estimate the relative contribution of nanoplankton (defined by fine POM in the present study) and microplankton (coarse POM) to the copepod diets, as well as copepods to the carnivore diets (Parnell et al., 2010(Parnell et al., , 2013. The model assumes that each isotopic ratio of consumers follows the pattern of a Gaussian distribution with an unknown mean and standard deviation. The structure of mean values of consumers is a weighted combination of the food sources' isotopic values. The 15 weights make up dietary proportions (given by a Dirichlet prior distribution). The standard deviation is divided up between the uncertainty around the fractionation corrections and the natural variability between all individuals within a defined group (Parnell et al., 2010(Parnell et al., , 2013. Because the values of consumers calculated from bulk copepod samples using the previous multiple linear regression models were only means and standard errors, we generated a vector consisting of 250 numbers for each group by a random normal distribution function. We then used the default iteration numbers (iterations = 500,000, burin = 50,000) 20 provided by the SIAR package to perform our analysis. Fractionation factors used in the model estimation were estimated by difference of TLs multiplying with 3‰ per TL for  15 N and with 0.5‰ per TL for  13 C, respectively. TLs were calculated from the  15 N difference between consumer and source as follows (Post, 2002): ( = 1 + ( 15 − 15 )/3).
Concentrations of isotope per mass among different diets (nanoplankton, microplankton, and major copepod genera) were not considered in this study. Model fitting was done via a Markov Chain Monte Carlo (MCMC) protocol that produces simulations 25 of plausible values of the dietary proportions of each source. More details on model simulation can be found elsewhere (Parnell et al., 2010(Parnell et al., , 2013. All statistical analyses were performed using R 3.4.0 software (https://cran.r-project.org/bin/windows/base/). Regression analyses of copepod isotopic values were performed by generalized additive models (GAMs) using the mgcv library (Wood and Wood, 2015). Data were smoothed by cubic regression splines and fitted by the family of Gaussian. One-way analysis of 30 variance (ANOVA) was adopted to test seasonal differences in environmental factors and copepod abundances, and Student's t-tests were used to test for significant differences in mean  13 C and  15 N values between nano-and microplankton. Before applying ANOVA and t-tests, the data were tested for normality of distribution and equal variance; significance was assumed at P = 0.05. For the coefficients of variation, we calculated them by dividing the standard error with mean value.

Environmental variability and zooplankton abundances
Environmental factors including temperature, salinity, Chl a levels, copepod abundance, dominant species, and percentages of 5 total copepods are shown in Table 1. Water temperature was significantly higher in summer and lower in winter (ANOVA, P < 0.001). Spatial variability of salinity was significant, with extremely low values at stations 1 and 2 (the river mouth) and then the values gradually increased to station 5 (the middle of the bay). Chl a concentrations ranged from 0.1 to 6.8μg l −1 and they were significantly higher in spring and summer than in winter and autumn (ANOVA, P < 0.01). The highest Chl a concentration occurred in the middle of the bay during the spring, while the lowest concentration was found at the river station 10 during the autumn. Seasonal variability of copepod abundance was significant (ANOVA, P < 0.01), with higher abundances in winter when temperatures and Chl a concentrations were low.
Detailed abundance composition of copepods is shown in Table S2. Seasonal and spatial variations of dominant species (>10% of total abundance) of copepods were apparent (Table 1). The marine calanoid Acartia dominated at the river mouth to the middle part of the bay, while Paracalanus dominated at the mouth of the bay during winter. Acartia also dominated at the 15 most highly saline stations in summer, except for station 7, where the community was dominated by Labidocera rotunda. A brackish water-preferring calanoid species, Pseudodiaptomus, dominated stations 1 and 2 at the river mouth in spring and another brackish calanoid species, Sinocalanus, dominated station 1 in autumn. At the river-mouth stations in summer, copepods were unexpectedly dominated by the marine calanoid species Tortanus dextrilobatus. The cyclopoid species Corycaeus affinis mainly dominated the most highly saline stations in spring and autumn. 20

Variability of plankton  13 C and  15 N values
The  13 C values of size-fractionated plankton (< 20 and 20-200 m) and mixed copepod samples showed distinct spatial variations in each season ( Fig. 2A-C). The  13 C values of nanoplankton (< 20 m POM) ranged from −27.6 to −19.4‰ with a mean of −22.7‰ (Fig. 2A). The lowest  13 C value of nanoplankton was found at station 1 (the upper stream station of Seomjin River) in spring and the highest at station 9 (the mouth of the estuary) in summer. The  13 C values of microplankton 25  m POM) ranged from −26.3 to −17.8‰ with a mean of −20.8‰ (Fig. 2B), being significantly higher than those of nanoplankton (paired t-test, t = 7.6, P < 0.001). Its lowest  13 C value was found at station 2 in spring and the highest at station 8 in winter. Overall, similar to nanoplankton, the microplankton  13 C values were more negative at the river portion (stations 1-3) and less negative at the mouth of the estuary (stations 7-9). The  13 C values of mixed copepods ranged from −25.9 to −16.4‰, with a mean of −20.1‰ (Fig. 2C). The lowest at station 1 in autumn and the highest at station 8 in summer. The 9 spatial variability of copepod  13 C values followed the pattern of POM  13 C values. However, the copepod  13 C values were significantly higher than those of nanoplankton (paired t-test, t = 8.6, P < 0.001) and microplankton (t = 3.1, P = 0.004). Their  13 C values were higher in summer and winter than in spring and autumn. At stations 1-3, river input lowered the  13 C values of nanoplankton during the wet season (spring to summer). At stations 4-9, significantly lower  13 C values were observed in autumn than in other seasons (ANOVA, F = 13.4, P < 0.001). For copepods, the autumn values were significantly lower than 5 those in other seasons (ANOVA, F = 5.9, P = 0.004).
The  15 N values exhibited wider fluctuations than  13 C values (coefficients of variation = 29.3% vs. 9.0%, 21.5% vs. 11.8%, and 18.8% vs. 13.1% for nanoplankton, microplankton, and copepods, respectively). The  15 N values of nanoplankton ranged from 3.2‰ (station 4 in summer) to 8.8‰ (station 1 in winter) with a mean of 5.6‰ (Fig. 2D). There were distinct patterns in the three locations of the bay. The  15 N values tended to decline with distance from the river mouth, then increased in the 10 middle of the bay, and decreased again toward the mouth of the estuary. The nanoplankton  15 N values were higher in winter than in other seasons (paired t-test, t = 5.4, P = 0.001 for spring; t = 3.0, P = 0.017 for summer; t = 4.1, P = 0.004 for autumn).
The  15 N values of the three major groups were all higher than the basal food resource (nanoplankton), and relatively higher than microplankton with some overlap of error bar. The values of harpacticoids isolated from the winter (stations 2-9) and 10 spring samples (stations 1) were measured directly. The mean  15 N values of harpacticoids (6.9 ± 0.6‰) were lower than those of other copepods and microplankton, but relatively higher  13 C value (−16.7 ± 1.4‰).
The trophic level in this study was defined as trophic position relative to nanoplankton, which was considered as the trophic 25 baseline. Considering the TL of nanoplankton is 1, the TL value of microplankton was calculated to be 0.7 times higher than that of nanoplankton (Fig. 4). As a whole assemblage balanced from different feeding behaviors, as indicated by the standard errors, copepods occupied a 1.2 level higher TL than that of nanoplankton, indicating herbivory (here, herbivory means a trophic level of 2) on nanoplankton with slight omnivory (TL = 2-3) on other dietary sources. The TLs of three major calanoid groups (marine calanoids and brackish calanoids) were similar to the bulk copepod assemblage with mean levels slightly higher 30 than 2. The mean TL value of cyclopoids had apparently higher trophic level than calanoids. In contrast, the mean TL of harpacticoids was very low, reflected by their low  15 N values.
Among calanoids, the mean TLs of Eurytemora pacifica (1 ± 0.5) and Pseudodiaptomus (1 ± 0.5) is indicative of their herbivorous and/or detritivorous characteristics. The mean TLs of Calanus sinicus (2.0 ± 0.6) and Centropages (1.6 ± 0.6) indicated they were primarily herbivorous; and those of Acartia and Paracalanus were slightly higher than 2 indicating they were primarily omnivorous feeding on both nanoplankton and microplankton. The levels of Sinocalanus (3.0 ± 0.7), Tortanus and Labidocera were higher than 3. 5 Based on TLs, the mean 15 N enrichments of the copepod assemblage were estimated to be 3.4‰ and 1.7‰ for nanoplankton and microplankton, respectively (Fig. 5A, B). The enrichment values on nanoplankton for both marine and brackish calanoids, as well as for the genera of Acartia were close to the average value of the copepod assemblage. The trophic enrichment of Calanus sinicus on nanoplankton (3.0 ± 1.7‰) was slightly lower than averaged marine calanoids and total copepods, while the enrichment on microplankton (1.7 ± 1.5‰) was similar. Centropages (1.9 ± 1.8‰) had much lower enrichment on 10 nanoplankton compared to total copepods and marine calanoids. Four high TL genera Tortanus, Labidocera, Sinocalanus and Corycaeus had high enrichments > 6 on nanoplankton and > 3 on microplankton.
The 13 C enrichments for total copepods were on average 0.6‰ and 0.3‰ when feeding on nanoplankton and microplankton, respectively (Fig. 5C, D). Patterns were same with enrichments on 15 N.The four high TL genera Tortanus, Labidocera, Sinocalanus and Corycaeus had high enrichments > 1.2 on nanoplankton and > 0.6 on microplankton. The enrichments of four 15 herbivorous/omnivorous species increased from Centropages, Calanus, and Acartia to Paracalanus. In contrast, the brackish calanoid genera Eurytemora and Pseudodiaptomus had extremely low enrichments of both 15 N and 13 C.

Contribution of size-fractionated POM to copepod diets
The Bayesian mixing model calculations showed that the contributions of different sizes of POM to copepod diets varied significantly with season (Student's t-test, P < 0.001 for all cases; Fig. 7). Size-selective feeding phenomena were particularly 20 apparent in winter (Fig. S2A) and summer (Fig. S2C). Mean contributions of microplankton accounted for about two-thirds of their assimilated diets at all stations in winter and summer, and were almost equal to that of nanoplankton in spring and autumn ( Fig. S2B, D). The proportions of the two size fractions of POM averaged from all four seasons contributing to copepod diets were also distinctly different at stations 2-5 and station 9 (Fig. S2). The mean contributions of microplankton to the copepod diets increased gradually from the river mouth up to a peak (0.81 ± 0.11) at the middle part of the bay. Then, the proportion 25 declined gradually to a trough (0.31 ± 0.18) at the deep-bay channel. The proportion then rebounded to a high level again at the bay mouth station.

Variability of  13 C and  15 N values of plankton with time and space
We found that seasonal variations and the spatial heterogeneity of copepod  13 C and  15 N values in Gwangyang Bay followed those of nanoplankton (POM < 20 m) and microplankton (POM > 20 m) (Figs 2 and 3). Based on the results of regression analyses, we found that the variability of copepod isotopic values was influenced by salinity (spatial variations), temperature 20 (temporal variations), and isotopic values of food sources (both spatial and temporal variations; Fig. 3). In general, spatial variations were much more pronounced because of the effect of river input and thereby riverine carbon in different salinity regimes. More negative values of three plankton groups (nanoplankton, microplankton, and copepods) were measured near the river mouth, and then the values increased progressively to the mouth of the estuary, indicating an apparently decreasing effect of river runoff and thus the uptake of carbon derived from river-borne terrestrial organic matter. These results are 25 consistent with other studies in estuarine environments (Cifuentes et al., 1988;Matson and Brinson, 1990;Thornton and McManus, 1994;Deegan and Garritt,1997;Fry, 2002). Such spatial distribution patterns have also been found for other primary producers such as seagrasses (reviewed by Hemminga and Mateo, 1996), macroalgae (Lee, 2000), as well as benthic microalgae (Kang et al., 2003), and the pattern will further propagate to consumers such as fish (Melville and Connolly, 2003;Herzka, 2005), oysters (Fry, 2002), mollusks (Antonio et al., 2010), and other benthic macro-invertebrates (Choy et al., 2008). 30 Seasonal successions of  13 C values were also apparent, probably because of high river input in spring, elevated productivity in summer, and species successions in autumn. When the wet season started in spring and phytoplankton started to bloom, river input lowered the  13 C values. The values increased again in summer because of the persistence of phytoplankton bloom (a low fractionation effect because of source limitation) and elevated productivity in summer. Although both river discharge and input of light carbon were low in autumn, the observed  13 C values were low. This was probably because of the lack of a 5 heavy carbon pool (a dissolved inorganic carbon pool in which the carbon was primarily composed of heavy 13 C) due to microbial respiration and species succession (Rau et al., 1990). During the post-bloom period in autumn, phytoplankton show low productivity and Chl a concentrations, with low abundance of diatoms but a dominance of flagellates (Baek et al., 2015).
Flagellates are known to have more negative  13 C values than those of diatoms arising from different fractionation effects (Gearing et al., 1984;Cifuentes et al., 1988;Rolff, 2000). 10 The  15 N variability of three major plankton groups was relatively complex spatially. Seasonal pattern of nanoplankton  15 N values, decreasing the values in summer to autumn, can well explained by the relationship between the  15 N values and temperature, indicated by significant GAM analysis. Spatial trends in the nanoplankton  15 N values can be explained by three distinct distribution patterns. The first pattern found in the river mouth area exemplifies a declining trend expected by mixing of freshwater planktonic materials, which grew up in water with high levels of dissolved inorganic nitrogen. The second pattern 15 found in the middle part of the bay, in which nitrate inputs were low while the concentrations of ammonia increased (Kwon et al., 2004; own data not shown), characterizes an increase in  15 N values in association with high Chl a concentrations.
Fractionations by autotrophic assimilation and bacterial utilization were the most likely source of the 15 N-enriched ammonia in nutrient pools of the middle of the bay (Cifuentes et al., 1988). The elevated POM  15 N values in the middle of the bay may be explained by 15 N-enriched ammonia remaining after algal uptake in the river mouth channel (Sato et al., 2006). The input 20 of sewage-derived 15 N-enriched ammonia (domestic sewage and livestock waste) could contribute substantially and swamp the other subtle processes to increase δ 15 N values of nanoplankton. The third distribution pattern represents declining  15 N values toward the offshore bay mouth in association with a reduction in the supply of 15 N-enriched nutrients from terrestrial sewage. Furthermore, the fractionation effect of phytoplankton will be reduced when nutrients substantially decreased and phytoplankton would take up nitrogen with little fractionation and stored relatively light nitrogen isotope (Cifuentes et al., 25 1988;Fogel and Cifuentes, 1993;Granger et al., 2004). As indicated by GAM regression analyses between the distribution of nanoplankton  15 N and environmental factors, significant increases in the nanoplankton  15 N values depend on ammonia and Chl a, further supporting our explanation. The microplankton  15 N values were stepwise elevated by environmental factors including temperature, salinity, ammonia and nitrate, indicated by significant results of GAM analysis (Table 2). Similar to nanoplankton, regression analysis also showed that the microplankton  15 N values increased significantly with increasing Chl 30 a concentrations (Table 2). One possible mechanism for this pattern is that higher phytoplankton abundance will result in a microplankton, indicative of the role of diatoms (preferring nitrate) in controlling the variation in microplankton  15 N values, whereas nanoflagellates (preferring ammonia) probably controlled the variation in nanoplankton  15 N values.
As indicated by GAM analysis, the seasonality of copepod 15 N values was primarily enhanced by temperature, which probably caused an elevated fractionation effect during the rapid assimilation of copepods. The lower explained deviance of  15 N indicated reflects the food-web processes that affect  15 N disproportionally and were not included in the GAM analysis. The 5 regression relationship between larger plankton and copepods was significant, whereas the patterns were somewhat decoupled, as they were primarily observed in spring and autumn. This kind of decoupling has also been reported in the open ocean (Montoya et al., 2002), where the transfer of nitrogen from primary producers to zooplankton is weak. A time lag in zooplankton development might cause the mismatch of zooplankton to 15 N-enriched POM at the initial stage of nutrient supplies. Indeed, here we found that the high  15 N values of copepods were primarily observed in summer, while the 10 corresponding  15 N values of POM started to increase from the winter, and phytoplankton blooming occurred in the spring.

Trophodynamics and trophic enrichments of copepods
The variability of copepod isotopic values in Gwangyang Bay suggests that the TLs of the copepod assemblage were highly dynamic. Because of different feeding behaviors and fractionation effects of copepods, the variability of the average Rhincalanus had the lowest values (Schmidt et al., 2003). A mesocosm study found that the  15 N values were increasingly higher in the order Temora < Pseudocalanus < Centropages, suggesting an increase of carnivory in the same manner (Sommer 20 et al., 2005). The trophic positions of primary consumers (Oithona and Neocalanus) and secondary consumers (Pleuromamma and Euchaeta) in the North Pacific Subtropical Gyre are estimated to be 2.1 and 2.9, respectively (Hannides et al., 2009). Furthermore, Kürten et al. (2011) reported that the relative trophic positions of zooplankton in the North Sea were high when the assemblage was mainly composed of Sagitta and Calanus, but low when the assemblage was dominated by Pseudocalanus and zoea larvae. 25 Our study has demonstrated the trophodynamics of estuarine copepods using multiple linear mixing model analysis based on the values of bulk samples and percentages in total biomass, by which the results of estimated  13 C and  15 N values were both significant (P < 0.01). The estimated  13 C values varied greatly up to 10.9‰ among groups (from the lowest for cyclopoids to the highest for harpacticoids) and 14.2‰ among genera (from the lowest for Paracalanus to the highest for Centropages), respectively (Fig. 4). The estimated  15 N values also varied somewhat up to 4.9‰ among groups, indicating one TL difference 30 among groups if we consider 3-4‰ trophic enrichment of  15 N between two adjacent TLs (Post, 2002;Michener and Kaufman, 2007). The trophic enrichments calculated by the differences from the basal food sources indicate that the overall enrichments of copepods were around 3.4‰ and 1.7‰ from nanoplankton and microplankton, respectively (Fig. 6). Our estimation shows that both marine calanoids and brackish water calanoids overall occupy a similar trophic niche (i.e., similar mean  15 N values), but have contrasting food sources (i.e., different  13 C values; Fig. 4A). The isotopic values of calanoid copepods indicate a mixture of different genera including both high and low  15 N values. For example, marine calanoids were mixed by high  15 N genera (Tortanus and Labidocera) and low  15 N genera (Acartia, Calanus, Centropages) and moderate  15 N genera 5 (Paracalanus and Acartia). Brackish types were mixed by high  15 N genera (Sinocalanus) and low  15 N genera (Eurytemora and Pseudodiaptomus). Consequently, the two major calanoids groups (marine and brackish water types) as well as the mixture of all copepod groups were estimated to be on average one TL higher than the nanoplankton base (Fig. 4). However, we do not necessarily conclude that they are herbivores because the nanoplankton studied here represent POM with a size range of 2-20 m, which may include ciliates, heterotrophic nanoflagellates, and heterotrophic dinoflagellates. Instead, all assemblages 10 mentioned above might be omnivorous with varying levels of relative trophic positions depending on dominant species.
Among calanoids, brackish water species had significantly lower  13 C values than marine species, indicative of an apparent effect of riverine carbon sources on brackish species through the food web (Fig. 3). These results were consistent with the distribution pattern that most brackish species occurred or dominant in the upper part of Gwangyang Bay. Sinocalanus tenellus had contrasting TL value to those of Pseudodiaptomus and Eurytemora (Fig. 4B), indicating a mixture of brackish water 15 calanoids being close to omnivory with a broad feeding size spectrum. However, as Sinocalanus and Pseudodiaptomus s were dominant in different seasons (Table 1), the TLs of the copepod assemblage at a specific condition will become relatively carnivorous (Sinocalanus dominating) or omnivorous (Pseudodiaptomus and Paracalanus parvus dominating). P. parvus is an important small species (body length 1mm) that is widely distributed in coastal and estuarine waters worldwide (Turner, 2004). Our results showed that, similar to other brackish water calanoids, this species was greatly influenced by 13 C-depleted 20 dietary sources, dominating both brackish stations in autumn and saline stations in winter (Table 1). This result indicates that Paracalanus was well adapted to fluctuating estuarine environments by feeding on prey originating from freshwater or prey that depends on riverine carbon sources. Acartia, one of the commonest genus in Gwangyang Bay throughout the year (Table   1), included both marine types (A. hudsonica and A. Omorii) and brackish types (A. ohtsukai, and A. erythraea) based on literature (Table S1). Marine types of Acartia primarily occurred or dominated during the winter and spring, while brackish 25 type dominated during the summer. Overall the genus of Acartia had higher  13 C values than those of Paracalanus, suggesting their different food sources. Acartia is also a widely distributed genus, with a switching feeding behaviour in response to the status of food composition (Kiøboe et al., 1996;Rollwagen-Bollens and Penry, 2003;Chen et al., 2013). The isotopic values of Acartia were similar to those of the assemblage of marine calanoids, indicating that this genus is omnivorous, as typical of marine calanoids. The marine calanoids species C. sinicus had a close trophic niche to total marine calanoids based on similar 30  15 N value and trophic level (Fig. 4), 1 TL higher than nanoplankton suggesting that this species is a typical marine herbivore relative to nanoplankton. Conversely, two other marine calanoid genera, Tortanus (T. dextrilobatus and T. forcipatus) and Labidocera (L. euchaeta and L. rotunda), were primarily carnivorous as indicated by their  15 N values (Figs. 3B and 4). These estimated results are consistent with the former experimental tests and field investigations (Ambler and Frost, 1974;Landry, 1978;Conley and Turner, 1985;Hooff and Bollens, 2004). Cyclopoids (primarily Corycaeus affinis) dominated copepod assemblages in spring and autumn at the middle part and deep-bay channel of the bay (Table 1). Our data reveal that Corycaeus was primarily carnivorous, being 2 TLs higher than nanoplankton (Fig. 4) and prefers feeding on 13 C-depleted dietary sources (Fig. 3). This result is consistent with previous reports that the Corycaeus genus is carnivorous (Gophen and Harris, 1981;5 Landry et al., 1985;Turner, 1984).
Isotopic values of microplankton indicate that they are roughly a half TL value higher than nanoplankton. Considering that the sizes of most ciliates and heterotrophic dinoflagellates primarily fall within this size spectrum (20-200 m), this result suggests an omnivorous trend among the mixed microplankton groups. Similarly, although measured only in winter samples, the benthic copepod group harpacticoids, represented by the species Euterpina sp., also differed from calanoids and cyclopoids with low 10  15 N values. The TL of harpacticoids estimated from this approach was somewhat misleading because of their unexpectedly low  15 N values, which probably reflect feeding on detritus or dead organisms that are depleted in 15 N (Sautour and Castel, 1993).

Selective feeding of copepods
Feeding preferences of different groups or genera on two size fractions of POM are of particular importance to explain seasonal 15 and spatial patterns of community trophic niches, and in turn will predict the impacts of the grazer community on lower TLs including phytoplankton and microzooplankton. Because not all possible food sources, such as bacteria, picoplankton, fecal pellets, and dead detritus, were investigated, our Bayesian mixing model calculations might have led to some biased results. Nevertheless, the model results might provide an estimation on what size fractions of dietary sources the grazer community ingest and assimilate. In general, our results highlight that the copepod assemblages have size-selective feeding behaviors, and 20 that these vary with season ( Fig. S2) and space (Fig. S3). The feeding selectivity of the bulk copepod assemblage was a balance of ingestion among different groups. The whole copepod assemblage assimilated two-thirds of its food requirement from microplankton in winter and summer, but they fed nearly equally on both size fractions of POM in spring and autumn (Fig.   S2). Our results suggest that groups that preferred large-sized prey played a more important role in total assemblage of copepods in winter and summer, during which the assemblage were primarily dominated by marine calanoids (Table 1). On 25 the other hand, in the spring and autumn when the assemblage was primarily dominated by carnivorous species (Corycaeus, Tortanus) and dominated by both Paracalanus (our result suggested this genus was an omnivorous species and the sizeselectivity was less pronounced preferring small-sized particle) and Acartia (preferring large-sized particle), the assemblage overall didn't show an apparent size-selectivity. Based on the model results for major copepod groups and genera, marine calanoids preferred feeding on large particles (Figs. 30 6A and 7A) Such a size-selective feeding preference has been widely reported in many field investigations (Liu et al., 2005a;Jang et al., 2010;Chen et al., 2017). On the contrary, brackish calanoids as a group preferred feeding on small-sized plankton.
These results were not well reported in literature. Among calanoids, Acartia and Paracalanus both contain marine species and brackish species. Our model results showed that Acartia preferred feeding on large-sized plankton, while Paracalanus was contrastingly different. Acartia, well studied in literature, are reported to prefer feeding on phytoplankton larger than 20 m in the coastal water adjacent to the present study area (Jang et al., 2010). Besides, due to nutritional requirements, Acartia were reported that they preferred predating microzooplankton such as ciliates and heterotrophic dinoflagellates (Wiadnyana and Rassoulzadegan, 1989;Chen et al., 2013). On the other hand, another dominant calanoid genus, Paracalanus, preferred feeding 5 on nanoplankton (Fig. 7E), suggesting that the feeding efficiency, with which the grazers' feeding appendage to retain particles, differs between these two genera. Because of such different feeding preferences, they dominated in different stations or seasons with differing food conditions (Table 1) and frequently co-occurred with little overlap of preferred food particle sizes. Besides, low  13 C value of P. parvus indicated that this wildly distributed species preferred feeding on particles flourished from terrestrial carbon. 10 Although C. sinicus was not dominant in abundance, this species contributes a great amount in biomass due to large body size (Table S1). Similar to Paracalanus, C. Sinicus preferred feeding on small-sized particles, while the difference is that the relatively high  13 C value of C. Sinicus indicated that this species preferred prey originated from offshore areas. The result is consistent with reports that C. Sinicus is a dominant species in Yellow Sea and East China Sea (the offshore shelf waters of this study estuary) (Uye, 2000;Liu et al., 2003). Another marine calanoid genus similar to Acartia, Centropages (C. 15 abdominalis and dorsispinatus) apparently preferred feeding on large-sized plankton (this study and Wiadnyana and Rassoulzadegan, 1989). For other two marine calanoid genera -Labidocera (L. rotunda and L. euchaeta) and Tortanus (T. dextrilobatus and T. forcipatus)feeding on two size fractions of POM did not occur, based on no contribution of the two size fractions of POM to their diets (Fig. 8A, C). Model results indicated that Labidocera and Tortanus were carnivorous genera, consistent with many reports (Landry, 1978;Mullin, 1979;Turner, 1984;Uye, 1994;Hooff and Bollens, 2004). Our 20 result also showed that Labidocera had an apparent selectivity in predating their prey, preferring marine calanoids (Acartia, Calanus and Centropages) and harpacticoids. This feeding selectivity was consistent with the distribution of this species. It never distributed in river stations (stations 1-3) (Table S2) and even dominated at station 7 during the summer (Table 1).
Differently, Tortanus didn't show selectivity among different prey (Fig. 8C) and thus this genus occurred frequently in different stations (Table S2). 25 The model results indicated that the size selectivity of brackish water calanoids such as Pseudodiaptomus and Eurytemora was also apparently for nanoplankton similar to Paracalanus (Fig.7 D, F). To our knowledge, the feeding habit of Pseudodiaptomus is unknown in the current literature, whereas some field studies suggest that estuarine Pseudodiaptomus flourishes by feeding on small phytoplankton cells (< 20 m) (Froneman, 2004), consistent with the present results. Incorporating the results of low  13 C values of Pseudodiaptomus and Eurytemora, we believe that these two genera were able to feed on those prey with small 30 sizes (2-20 μm) and low  13 C from originated from terrestrial sources, whereas low trophic position and trophic enrichments indicated that these two genera may ingest prey with extremely low  15 N compared to nanoplankton (e.g., detritus). The detritus with low δ 15 N may contribute to the balance the δ 15 N of Pseudodiaptomus and Eurytemora, and thus, they may act as detritivores besides of herbivores. Another brackish calanoid species -Sinocalanus tenellus were more diverse in feeding selectivity (Fig. 8B). It can act as a suspension feeder preferring small-sized plankton as well as a raptorial feeder preferring predating other brackish calanoids and cyclopoids. This species was even reported as a cannibalistic feeder in estuary (Hada and Uye, 1991).
In contrast to both marine calanoids and brackish water calanoids, the cyclopoid species Corycaeus affinis was a carnivorous 5 species preferring predating on two brackish species -Paracalanus and Pseudodiptomus (Fig. 8D). This species is a widely distributed in all the world's oceans (Turner, 2004) and also frequently dominate in copepod assemblages in spring and autumn in Gwangyang Bay (Table 1). Although the isotopic data for harpacticoid copepods were limited to only one season, we still obtained a clear feeding selectivity pattern based on the Bayesian mixing model (Fig. 7C). Harpacticoids in winter preferred microplankton to nanoplankton if benthic food sources were not considered. Their feeding selectivity contributed to the overall 10 feeding preference of total copepods in this season.

Conclusions
Here we have demonstrated the temporal and spatial variability of stable isotope ratios of copepods, which was determined by the isotopic values of two size fractions of POM, and strongly influenced by salinity (spatiality) and temperature (temporality).
Such characteristic is a key in understanding the biogeochemical cycles of carbon and nitrogen in Gwangyang Bay. We further 15 used a simple linear mixing model and a Bayesian mixing model to extrapolate from the information derived from the isotopic analysis of bulk copepod samples. The model results were robust and allowed the estimation of the relative TLs and trophic enrichment (fractionation effect) of different groups and dominant genera of copepods, as well as their diet compositions.
Temporal and spatial patterns of copepod isotopic traits were further explained by size selectivity on plankton size fractions, as well as the feeding preference of dominant species. Based on such relative trophic positions and feeding preference, we can 20 depict a simple energy flow of the Gwangyang Bay planktonic food web: from primary producers (nanoplankton) and a mixture of primary producers and herbivores (microplankton) through omnivores (represented by Calanus, Centropages, Acartia, and Paracalanus) and detrivores (represented by Pseudodiaptomus and Eurytemora) to carnivores (dominated by Corycaeu, Tortanus, Labidocera, and Sinocalanus) (see a simplified energy flow, Fig. 9).  2: Coefficients (F) and significance levels (p) of the effects of variable predictors on δ 15 N of nanoplankton (plankton < 20 µm) and microplankton (plankton > 20 µm) using a Generalized Additive Model test. P-value < 0.05 indicates significance. The symbol n.s. represents no significance.