Articles | Volume 18, issue 2
Research article
07 Mar 2022
Research article |  | 07 Mar 2022

Carbon and nitrogen dynamics in the coastal Sea of Japan inferred from 15 years of measurements of stable isotope ratios of Calanus sinicus

Ken-ichi Nakamura, Atsushi Nishimoto, Saori Yasui-Tamura, Yoichi Kogure, Misato Nakae, Naoki Iguchi, Haruyuki Morimoto, and Taketoshi Kodama

Both nitrogen and carbon dynamics have changed in the Sea of Japan. We hypothesized that the carbon and nitrogen stable isotope ratios (δ13C and δ15N) of the copepod Calanus sinicus could record changes in the coastal environment of the Sea of Japan. Consequently, these isotope ratios were monitored during the spring at four stations from 2006 to 2020 to identify the changes in carbon and nitrogen dynamics. The δ13C values ranged from 24.7 ‰ to 15.0 ‰ and decreased from the spring bloom (February–March) to the post-bloom (June–July) seasons. These variations were attributed to changes in the physiology of both C. sinicus and phytoplankton δ13C contents. The δ15N values range from 2.8 ‰ to 8.8 ‰, indicating that C. sinicus is a secondary producer; the tendency of the δ15N values to increase from the bloom to the post-bloom seasons was attributable to an increase in the δ15N of phytoplankton. A generalized linear model (GLM) approach indicated that >70 % of the variations in δ13C can be explained by sea surface temperature (SST), sea surface chlorophyll a concentration (SSC), carbon:nitrogen ratio of C. sinicus (C/N ratio), and geographic differences. The residuals of δ13C in the GLM decreased yearly (0.035 ‰ yr−1). The GLM for δ15N of C. sinicus indicated that δ15N varies with the stage or sex in addition to SST, SSC, C/N ratio, and geographic differences. The δ15N values of female C. sinicus and stage V copepodites were the lowest and highest, respectively. The residuals of δ15N in the GLM did not exhibit a significant interannual trend. These results suggest that the carbon isotope ratio in the secondary producer has linearly changed in the coastal Sea of Japan over the past 15 years. Moreover, the changes in carbon dynamics of this area have been recorded and observed to impact the marine ecosystem, while the nitrogen dynamics have not been recorded despite the increasing nitrogenous nutrient inputs in this sea.

1 Introduction

Coastal ocean ecosystems are ubiquitously important and have greatly changed as a result of human activities (Halpern et al., 2008; Doney, 2010). Stable isotope ratios of carbon and nitrogen have been employed to discern both the elemental cycles and environmental changes in marine ecosystems (Ohman et al., 2012; Lorrain et al., 2020; Ren et al., 2017), along with metrics of the trophic positions of organisms in an ecosystem (Aita et al., 2011). During the 21st century, the 13C:12C carbon isotopic ratio of tuna has rapidly been decreasing as a linear function with time, and the rate of decrease has been faster than expected owing to the Suess effect and the increase in anthropogenic carbon emissions (Gruber et al., 1999). Thus, it has been hypothesized that lower-trophic-level ecosystems have changed globally (Lorrain et al., 2020). In the marginal seas of East Asia, the 15N:14N nitrogen isotopic ratio of organic matter bound in coral skeletons has decreased with an increase in anthropogenic nitrogen deposition (Ren et al., 2017).

The Sea of Japan (Japan Sea) has been significantly impacted by human activities and global climate change (Chen et al., 2017). Monotonic changes have been detected in the chemical environment of coastal areas of the Sea of Japan, as well as in the global oceans (Ishizu et al., 2019; Ono, 2021; Kodama et al., 2016). In the surface waters of the Sea of Japan, pH and concentrations of phosphate and dissolved oxygen have been decreasing during the last few decades owing to global warming and the increase in atmospheric carbon dioxide (Ishizu et al., 2019; Ono, 2021; Kodama et al., 2016). The anthropogenic nitrogen inputs from the atmosphere have been increasing in the Sea of Japan (Kitayama et al., 2012), as well as in the global oceans (Duce et al., 2008). The Sea of Japan is regarded as an oceanic microcosm of the changing global ocean (Chen et al., 2017); thus, the impacts of changing inorganic chemical environments on biological production can be clearly detected compared to those in other oceans. Notably, the muscles of small pelagic fish in the Sea of Japan and the East China Sea from 1996 to 2019 showed that the annual mean 13C:12C and 15N:14N values monotonically decreased by 0.08 ‰ yr−1 and 0.05 ‰ yr−1, respectively (Ohshimo et al., 2021). This suggested a decrease in the stable isotope ratios of carbon and nitrogen in the food of small pelagic fish in the Sea of Japan. However, isotope values of both zooplankton and phytoplankton, which reflect the chemical environment in the Sea of Japan, have rarely been reported.

Calanus sinicus, the focus of this study, is a dominant copepod species in the coastal waters of the western North Pacific (Uye, 2000), including the Sea of Japan (Kodama et al., 2018a). Copepods of the genus Calanus are the major components of macrozooplankton in shelf and coastal ecosystems outside the tropics and are the major source of nutrition for pelagic fish (Uye, 2000). Maximum longevity of C. sinicus depends on the temperature but is usually less than 1 year (Uye, 1988); the development time from hatching to adulthood is ∼50 d at 10 C and 16 d at 20.2 C (Uye, 1988). The reproduction of C. sinicus occurs throughout the year but mainly at the end of spring or early summer (Huang et al., 1993; Zhang et al., 2005). The diapause phase is observed in the summer in the Yellow Sea (Pu et al., 2004b) but not in the Seto Island Sea of Japan (Huang et al., 1993). Calanus stores a large amount of lipid reserves during phytoplankton blooms (Gatten et al., 1980). The biology and ecology of C. sinicus in the Sea of Japan have been poorly studied, but the coastal areas of the Sea of Japan are spawning and nursery grounds for Japanese sardines (Sardinops melanosticta), Japanese anchovy (Engraulis japonicus), and Pacific bluefin tuna (Thunnus orientalis) (Nishida et al., 2020; Furuichi et al., 2020; Ohshimo et al., 2017); the larvae of these species prey on both the larval and adult stages of C. sinicus (Hirakawa et al., 1997; Hirakawa and Goto, 1996; Kodama et al., 2017a). Therefore, C. sinicus plays a vital role as a major component of the second trophic level in the Sea of Japan.

In this study, we hypothesize a decrease in carbon and nitrogen stable isotope ratios of C. sinicus in addition to that in small pelagic fish (Ohshimo et al., 2021), as the changing chemical environments in the Sea of Japan affected the marine ecosystem. Hence, we investigated the carbon and nitrogen stable isotope ratios of C. sinicus in the coastal area of Japan for 15 years in the 21st century and determined the controlling factors of stable isotope ratios using an empirical model approach. Subsequently, we calculated the long-term trends in the stable isotope ratios of C. sinicus to detect the impacts of the changing chemical environments in the Sea of Japan on the marine ecosystem.

2 Materials and methods

2.1 Onboard observations

Onboard observations were conducted from the pre- to post-bloom seasons from 2006–2020 during the R/V Mizuho-Maru, R/V Yoko-Maru (Japan Fisheries Research and Education Agency), and R/V Dai-Roku Kaiyo-Maru (Kaiyo Engineering Co., Ltd.) cruises in the territorial waters in the Sea of Japan, a marginal sea of the western North Pacific (Fig. 1). Four stations were chosen for the collection of stable isotope subsamples: Toyama Bay (TB), Iida Bay (IB), north of the Noto Peninsula (NN), and Wakasa Bay (WB). These four stations are among the 26 stations described in a previous study (Kodama et al., 2018a) that reported a clear west–east gradient of the zooplankton community structure with a temperature gradient in this area during May. Based on a previous study (Kodama et al., 2018a), the TB and IB stations were located in a cold-water area, and the NN and WB stations were located in a warm-water area. Additionally, the TB station was located near the mouths of two rivers (the Sho and Oyabe rivers; Fig. 1), and the WB station was in an area of restricted circulation. The cruises were conducted at four time intervals: (1) the end of February and/or the beginning of March (described as March), (2) the end of April, (3) the middle or end of May, and (4) the end of June or beginning of July (described as June) (Fig. 1c). The cruise in March corresponded to the early stage of the spring phytoplankton bloom; that in the April occurred during the late stage of the bloom, and those in May and June occurred during post-bloom conditions, according to Kodama et al. (2018b). Observations after 2015 were conducted only in April and May (Fig. 1c).

Figure 1Sampling locations and dates. (a) Location of the study area in the Japan Sea. (b) Distribution of sampling stations of Calanus sinicus (closed circles) and particulate organic matter (POM) for stable isotope analysis along the Japanese coast of the Japan Sea. (c) Months and years of samplings of C. sinicus at the four stations. The lines in (b) denote the downstream reaches of class A rivers. The maps were made with Natural Earth and the Geospatial Information Authority of Japan.

Copepods for stable isotope analyses were collected by towing a bongo net (500 µm mesh and 70 cm mouth diameter) obliquely at 0.5 s−1 from a depth of 75 m, or 10 m above the bottom, to the surface at each station. Collected samples were frozen and preserved at -20 C until sorting was conducted at an onshore laboratory. Some samples were stored for more than 10 years. Vertical profiles of temperature, salinity, and chlorophyll a concentration were measured using a conductivity–temperature–depth (CTD) sensor (Sea-Bird, SBE 9 plus, or SBE 19 plus) with an in vivo chlorophyll fluorescence sensor from a depth of 200 m to the surface. The temperature and salinity of the sensors were annually calibrated by the manufacturer. Surface seawater was sampled with a bucket to measure the sea surface temperature (SST), sea surface salinity (SSS), and sea surface chlorophyll a concentration (SSC). SST and SSS were measured with a calibrated mercury thermometer and salinometer (Autosal, Guildline Instruments), respectively. To measure the chlorophyll a concentrations, particles in 300 mL of water were collected on a glass fiber filter (GF/F, Whatman), and the pigments were extracted with N,N-dimethylformamide. Chlorophyll a concentrations were estimated based on the fluorescence of the extract, which was measured using a fluorometer (10-AU, Turner Designs) (Holm-Hansen et al., 1965). The chlorophyll fluorescence sensor was calibrated using discrete samples collected during each cruise. The nutrient concentrations at the surface were measured in cruises conducted after 2015. The procedures used for the nutrient analyses have been described by Kodama et al. (2015).

We collected and measured the stable isotope ratios of particulate organic matter (>0.7µm, POM) in April 2017 (at 0 m depth, n=6) and May 2019 (at 10 and 30 m depths, n=6) near the C. sinicus sampling sites (Fig. 1b). Particles in 1.5–4.6 L seawater were collected on the pre-combusted (460 C, 6 h) glass fiber filter (GF/F, Whatman), which were frozen and stored until the analyses in the onshore laboratory.

All data and program codes are obtained from the Digital Commons Data (Nakamura et al., 2021).

2.2 Stable isotope analyses

A total of 94 frozen Bongo net samples were thawed at room temperature. From each thawed sample, Calanus sinicus (Copepoda; Calanoida) adults or stage V copepodites were sorted as quickly as possible under a dissecting microscope to avoid their alteration or contamination. Calanus sinicus samples were soaked in a 3 % NaCl aqueous solution to clean their bodies. The 28 samples were divided into subsamples of stage V copepodite (C5), adult female (F), and adult male (M), and the remaining 66 samples were not divided. The physiological processes of C. sinicus differ according to their stage and sex (Zhou et al., 2016; Pu et al., 2004a, b). At least five individuals were necessary for measuring the stable isotope ratio in each subsample. When C. sinicus was present in one sample, several subsamples were prepared to evaluate the variability of the analysis. Therefore, a total of 267 subsamples were prepared to measure the stable isotope ratios. Each subsample was wrapped in a tin disk after drying at 60 C in a drying oven for 36–48 h. Table 1 presents the details regarding the number of subsamples. The stable isotope ratios of each subsample were measured.

Table 1Numbers of subsamples of Calanus sinicus for stable isotope analysis.

Download Print Version | Download XLSX

The carbon and nitrogen stable isotope ratios of the subsamples were measured using a stable isotope ratio mass spectrometer (IsoPrime100; Elementar) coupled with an elemental analyzer (vario MICRO cube, Elementar). The stable isotope ratios of carbon and nitrogen were calculated in per mill (‰) deviations from the corresponding standards using the following equation:

(1) δ 13 C bulk or δ 15 N bulk = [ ( R sample / R standard ) - 1 ] × 1000 ,

where R is the 13C/12C or 15N/14N ratio, respectively. The standards for carbon and nitrogen were L-alanine, and the references were Pee Dee Belemnite and atmospheric N2, respectively. The precision of the analysis was within 0.2 ‰ for both δ13C and δ15N. Based on the precision of the analyses of L-alanine, the δ13C and δ15N values for each subsample were rounded off to one decimal place. The amounts of carbon and nitrogen in each subsample were also recorded, and the carbon:nitrogen ratio (C/N ratio) of C. sinicus was calculated. However, we did not conduct any lipid extraction processes. Lipids affect δ13C; however, lipid corrections are possible using modeling approaches (Logan et al., 2008). Thus, we calculated the lipid-free δ13C (δ13Cex) based on the δ13Cbulk, C/N ratio, and the equation reported by Smyntek et al. (2007).

Before measuring the isotope ratios of POM, the filters were exposed to HCl fumes for more than 2 h to remove carbonate salts. The decarboxylated filters were dried in an oven (60 C, overnight), and the entire filter was wrapped in a tin cup to measure the carbon and nitrogen stable isotope ratios, as well as the C. sinicus.

Figure 2Variations of environmental parameters among months and stations. (a) Sea surface temperature (SST), (b) sea surface salinity (SSS), and (c) sea surface chlorophyll a concentration (SSC). Values at TB, IB, NN, and WB are shown. Box plots show the medians (horizontal lines within boxes), upper and lower quartiles (boxes), and quartile deviations (bars). Small points are raw data points. In (a), the lowest outlier at TB in March (SSS = 28) is not shown.


2.3 Statistical analyses

The relationships between δ13Cbulk and δ15Nbulk of C. sinicus and environmental parameters (SST, SSC, and SSS) were analyzed using general linear models (GLMs) with the linear link function using the R software (R Core Team, 2020). The structures of GLMs were based on the concepts articulated by Kodama et al. (2021). The δ13C and δ15N errors were assumed to be normal distributions, and explanatory variables were expressed as linear functions in the GLMs, with the exception of SSC, which was log-transformed. Two types of GLMs were applied as follows:

(2) δ 13 C bulk or δ 15 N bulk poly ( SST , 2 ) + poly ( log SSC , 2 ) + poly ( SSS , 2 ) + poly ( C / N , 2 ) + f ( stn ) + f ( stage ) ,

where stn and stage represent the target δ13Cbulk or δ15Nbulk, station (WB, NN, IB, and TB), and copepod stages (C5, F, M, and mixed), respectively. The arguments of the f functions are categorical variables used to simulate nonlinear relationships. The second argument of the polyfunctions (i.e., 2) indicates that the first one was incorporated into a quadratic expression in the model; namely, dome-shaped responses were included in the models. We added the “mixed” stage category for the subsamples undivided into stage or sex categories. We applied this GLM approach to δ13Cex; however, the results, except for the C/N ratio, were very similar to δ13Cbulk, and thus the results of δ13Cex are not shown. The significance of the parameters was checked by the analysis of variance (ANOVA). The explanatory variables in the final version of Eq. (2) were determined based on the values of generalized variance inflation factors (GVIFs) and Akaike information criteria (AICs). We required the GVIFs of the explanatory variables to be less than 10, and the model with the smallest AIC was accepted as the final model. For example, we did not apply f (month) and f (year) to the explanatory variables in Eq. (2) because they are collinear with other environmental parameters. The “ggpredict” function in the “ggeffects” package (Lüdecke, 2018) was used to visualize the effect of the explanatory variables on the δ13Cbulk or δ15Nbulk values. When categorical variables (i.e., stage and station) remained in the models with the least AIC, their values were fixed to calculate the predicted values.

We evaluated the monthly and interannual variations of the stable isotope ratios using the residuals (observation–prediction) of Eq. (2) (described as residual δ13C and δ15N). The trends of the residual δ13C and δ15N were expected to indicate the trend of δ13Cbulk and δ15Nbulk, excluding the interannual trends of environmental parameters and geographical effects.

3 Results

3.1 Environmental variables

The SST increased from March to June at every station (Fig. 2a). The SST at WB, the westernmost station, showed the highest SST among the four stations during every cruise. The SST ranges were 8.8–23.1, 9.5–21.9, 9.3–23.2, and 10.1–23.4 C at stations TB, IB, NN, and WB, respectively. The SSS was low at station TB (Fig. 2b). The ranges of SSS were 28.6–34.5, 33.1–34.3, 33.1–34.6, and 33.1–34.6 at stations TB, IB, NN, and WB, respectively. The SSS was the highest in May at the other three stations. Variations of SSC concentration at TB differed from those at the other stations. The SSC was higher at TB than at the other stations in every month. Additionally, the SSC at TB was higher in June than in May (Fig. 2c). The SSC values ranged from 0.02–7.45, 0.08–5.46, 0.05–9.07, and 0.08–6.60 µg L−1 at stations TB, IB, NN, and WB, respectively.

The mean ± standard deviation of the C/N ratio in C. sinicus was 5.5±2.6 g g−1. The C/N ratio in C. sinicus ranged from 3.4 to 19.2 g g−1. The C/N ratio was significantly different among the stages (ANOVA, p<0.001, degree of freedom = 3, F = 30.1); the C/N ratio was highest in C5 (mean ± sd: 7.9 ± 4.0 g g−1, n=61) and lowest in F (4.6±0.8 g g−1, n=78). The monthly variations were also significant (ANOVA, p<0.001, degree of freedom = 3, F = 15.09); the C/N ratios in all stages were the highest in April (Fig. 3a).

Figure 3Variations of (a) C/N ratio, (b) δ13Cbulk, (c) δ13Cex, and (b) δ15Nbulk values in C. sinicus among growth stages and sex. The growth stages and sex of C. sinicus were divided into copepodite V (C5), adult female (F), adult male (M), and mixed (undivided). Box plots show the medians (horizontal lines within boxes), upper and lower quartiles (boxes), and quartile deviations (bars). Small points are raw data points, and the colors and shapes of points indicate stations and sampling months, respectively.


3.2 Stable isotope ratios and their relationship between environmental parameters

The δ13C of POM ranged from 27.3 ‰ to 22.5 ‰ (mean ± standard deviation: -24.8±1.6 ‰) and from 24.9 ‰ to 23.0 ‰ (23.7 ‰ ± 0.7 ‰) in April 2017 and May 2019, respectively. The δ15N of POM ranged 3.1 ‰ to 5.1 ‰ (4.1 ‰ ± 0.7 ‰) and 2.5 ‰ to 4.3 ‰ (3.5 ‰ ± 0.7 ‰) in April of 2017 and May of 2019, respectively.

The means ± standard deviations of all of the C. sinicus δ13Cbulk, δ13Cex, and δ15Nbulk values were 20.6 ‰ ± 1.8 ‰ (24.7 ‰ to 15.0 ‰), 19.6 ‰ ± 1.9 ‰ (24.8 ‰ to 14.7 ‰), and 6.9 ‰±1.2 ‰ (2.8 ‰ to 8.8 ‰), respectively. The highest δ13Cbulk values were observed in March at all four stations (TB: 17.5 ‰ ± 2.1 ‰; IB: 18.7 ‰ ± 1.1 ‰; NN: 18.7 ‰ ± 0.9 ‰; and WB: 18.5 ‰ ± 1.1 ‰; Fig. S1 in the Supplement), and the lowest monthly mean δ15Nbulk was observed in March at all four stations (TB: 5.1 ‰ ± 0.8 ‰; IB: 5.1 ‰ ± 1.4 ‰; NN: 5.6 ‰ ± 1.3 ‰; and WB: 4.5 ‰ ± 1.0 ‰; Fig. S1). Comparing the δ13Cbulk values among the C. sinicus stages, those of C5 were significantly lower than those of F and M (Fig. 3b, Tukey's honestly significant difference – HSD – test, p<0.001), whereas no significant differences were observed in the δ13Cex values between the stages (Fig. 3c, ANOVA, p=0.186). Comparing the δ15Nbulk values among the C. sinicus stages, those of C5 were significantly higher than those of F and M (Fig. 3d, Tukey's HSD test, p<0.017), and those of F, M, and mixed stages were at the same level (Tukey's HSD test, p>0.8). The δ13Cbulk, δ13Cex, and δ15Nbulk values of C. sinicus in April 2017 were 21.6 ‰ ±  0.6 ‰, 18.2 ‰ ± 1.3 ‰, and 8.1 ‰ ± 0.6 ‰, respectively; those in May 2019 were -23.7±0.1 ‰, -22.2±0.1 ‰, and 6.7 ‰±0.4 ‰, respectively.

In the model with the least AIC, we used Eq. (2) to remove SSS from the full model for both δ13Cbulk and δ15Nbulk. Additionally, this stage was not selected in the least-AIC model of δ13Cbulk. Thus, the following equations represent the least-AIC models.


The coefficient of determination (r2) values of the least-AIC δ13C and δ15N models were 0.712 and 0.448, respectively. The effects of all variables were significant, except for the logSSC stage in δ15Nbulk (ANOVA). The dispersion values of both models were <1.

The least-AIC models produced convex graphs of δ13C as a function of SST or the C/N ratio (Fig. 4a, c). The minima of δ13Cbulk values occurred at an SST and a C/N ratio of approximately 20 C and 14, respectively (Fig. 4a, c). Inter-station comparisons indicated that the δ13Cbulk values were 0.6 ‰–1.1 ‰ higher in WB and TB compared to those in IB and NN (Fig. 4d).

Figure 4Responses of least-Akaike information criterion (AIC) generalized linear models (GLMs) based on least-squares mean (predicted) values. (a) δ13Cbulk versus sea surface temperature (SST), (b) δ13Cbulk versus sea surface chlorophyll a concentration (SSC), (c) δ13Cbulk versus carbon:nitrogen (C/N) ratio, (d) δ13Cbulk versus stations, (e) δ15Nbulk versus SST, (f) δ15Nbulk versus SSC, (g) δ15Nbulk versus C/N ratio, (h) δ15Nbulk versus station, and (i) δ15Nbulk versus stages. Black lines with gray shading (ac, eg) or closed circles with vertical bars (d, h and i) denote the predicted values with 95 % confidence limits based on the GLMs.


The responses of δ15Nbulk to SST mirror those of δ13Cbulk to SST (Fig. 4e). The least-AIC model produced a concave graph of δ15Nbulk as a function of logSSC (Fig. 4f). The δ15Nbulk maxima occurred at SST and SSC values of approximately 18 C and 1 µg L−1, respectively (Fig. 4e, f). Inter-station comparisons indicated that δ15Nbulk values were lowest in WB and highest in NN (Fig. 4g). The comparison among the stages indicated that δ15Nbulk values were lowest in F and highest in C5 and the mixed subsamples (Fig. 4h).

3.3 Temporal variations of residuals

The residual δ13C (residuals of δ13Cbulk based on the GLM described by Eq. 3) showed significant interannual variations (ANOVA, p<0.001, Fig. 5a) but no significant monthly variations (ANOVA, p=0.09, Fig. 5b). Additionally, the residual δ13C showed a significant negative relationship between the years (t test, p=0.0035, Fig. 5a). The coefficient (± standard error) of the relationship between residual δ13C and year was 0.0365 ‰ yr−1±0.0124 ‰ yr−1. When calculating the relationships between the residual δ13C and year at every station, significant negative relationships were only observed at the NN station (t test, p<0.001). High annual mean residual δ13C values were observed in 2006 (mean ± standard deviation : 0.51 ‰±0.89 ‰), 2007 (0.72 ‰ ± 0.99 ‰), 2015 (0.44 ‰ ± 1.01 ‰), and 2020 (0.55 ‰ ± 0.11 ‰) (Fig. 5a).

Interannual variations in residual δ15N (residuals of δ15Nbulk based on the GLM described by Eq. 4) were also significant (ANOVA, p<0.001; Fig. 5c), but the linear trend was not significant (t test, p=0.053). The monthly variations in residual δ15N were significant (ANOVA, p<0.001) and increased significantly from March to April (Fig. 5d). Low annual mean residual δ15N values were observed in 2012 (0.59 ‰ ± 0.63 ‰), 2014 (0.80 ‰ ± 0.81 ‰), and 2020 (0.97 ‰ ± 0.73 ‰) (Fig. 5a).

4 Discussion

We hypothesize that the δ13C and δ15N values of C. sinicus decreased over time. Our GLM approach indicated that the δ13C, and not the δ15N, decreased yearly. To confirm that this trend was attributed to the chemical environment, we discussed the trophic position of C. sinicus in the Sea of Japan and the effects of lipids. Finally, we compared our results with those of other studies to identify the characteristics of the chemical environment in the Sea of Japan.

Figure 5(a) Interannual and (b) monthly variations of residual δ13C (residuals of δ13C in the GLM) values and (c) interannual and (b) monthly variations of residual δ15N values. Box plots show the medians (horizontal lines within boxes), upper and lower quartiles (boxes), and quartile deviations (bars). Small points are raw data points, and the colors of points indicate stations. The black lines and gray shading (a, c) indicate linear regression lines and 95 % confidence intervals, respectively. The coefficient of regression line was significantly different from 0 between residual δ13C and year (solid line) but not between residual δ15N and year (dotted line).


In this study, we assumed that C. sinicus is the “secondary producer” in the Sea of Japan because C. sinicus is known to prey on heterotrophic plankton in addition to phytoplankton (Hirai et al., 2018; Yi et al., 2017; Uye and Yamamoto, 1995). The mean δ15Nbulk values of C. sinicus (6.9 ‰ ± 1.2 ‰) are consistent with this scenario; their δ15N values are intermediate between those of POM in coastal areas of the Sea of Japan (around 2 ‰–6 ‰; Kogure, 2004; Antonio et al., 2012) and predatory fishes such as Japanese anchovy (8.9 ‰–11.7 ‰; Tanaka et al., 2008) and Japanese sardines (9.4 ‰ ± 0.7 ‰; Ohshimo et al., 2019). The differences in δ15N between POM and C. sinicus were 4.0 ‰ and 3.2 ‰ in April 2017 and May 2019, respectively. These differences are within the range of the nitrogen discrimination factor of 3.0 ‰ ± 1.0 ‰ (Aita et al., 2011). We found that δ15Nbulk values were different among the stages (CF, F, and M), suggesting that the prey or nitrogen metabolism differed among the stages; however, the differences among the stages were small (≤0.5 ‰) and not significant in the GLM analysis. Therefore, the case studies in April 2017 and May 2019 suggested that the δ15Nbulk of C. sinicus reflects the δ15N values of POM.

In contrast to δ15N, the relationship between δ13Cbulk of C. sinicus and δ13C of POM differed between April 2017 and May 2019. The δ13Cbulk of C. sinicus was 3.0 ‰ higher than that of POM in April 2017, but the difference was 0 ‰ in May 2019. C. sinicus can store oil in the sac (Zhou et al., 2016). The δ13C values of copepods decrease with increasing fatty acid content (Smyntek et al., 2007). However, the relationship between δ13Cex of C. sinicus and δ13C of POM was different between April 2017 (6.6 ‰) and May 2019 (1.5 ‰), similar to the trends of δ13Cbulk in the two periods. The δ13C of POM was not very different between April 2017 and May 2019, and the spatial variations were not as large. The in situ δ13C of POM, however, is not immediately reflected in the δ13C of C. sinicus; the turnover of carbon in copepod tissues must also be considered. The turnover of Calanus finmarchicus has been reported to be 5–10 d (Mayor et al., 2011), and the development time of C. sinicus is similar to that of C. finmarchicus (Uye, 1988). Considering the turnover time of C. sinicus and the spatiotemporal heterogeneity of δ13C of POM, neither δ13Cbulk nor δ13Cex of C. sinicus reflected the δ13C of POM at the same station. Aita et al. (2011) reported that the relationship between δ15N and δ13C in marine plankton ecosystems differed between seasons in the Oyashio area, western North Pacific, and the slopes of linear regression between δ15N and δ13C varied between 0.61 (May) and 1.39 (July). Therefore, variations in δ13Cbulk and δ13Cex of C. sinicus not only reflected the variations in δ13C of POM, but also highlighted the carbon discrimination factor and metabolism of C. sinicus. Neither δ13Cbulk nor δ13Cex reflected the δ13C of POM; however, δ13Cbulk was well explained (r2=0.712) with the environmental parameters of the GLM approach using Eq. (3). This is not questionable because the quantities of temperature and food (i.e., chlorophyll a concentration) are known to control the discrimination and metabolism of C. sinicus (Jager et al., 2015).

The C/N ratio was different among the growth stages and sexes, thereby suggesting that lipid contents were different among them. High amounts of lipids decrease the δ13Cbulk of zooplankton (Smyntek et al., 2007). Modeling approaches are usually applied to remove lipid effects (Logan et al., 2008). The GLM results in our study indicated that considering the environmental parameters, the δ13Cbulk of C. sinicus was not different among the growth stages of copepods. In our study, the model did not completely remove the effects of lipids, but the results of GLM suggested that the effect of lipids on residual δ13C was negligible and was treated the same as lipid-extracted δ13C.

A significant linear decreasing trend in the residuals of δ13Cbulk (0.0365 ‰ yr−1± 0.0124 ‰ yr−1) was observed in our study. The decreasing trend of δ13Cbulk was slower than that of δ13C values of tuna muscle (0.12 ‰ yr−1; Lorrain et al., 2020) but was the same as that of δ13C in small pelagic fish in the Sea of Japan and East China Sea: from −0.05 % yr−1 to −0.08 % yr−1 (Ohshimo et al., 2021). Additionally, our decreasing trend was comparable to the Suess effect, which was reported to be 0.025 ‰ yr−1 (Gruber et al., 1999). This suggests that the decreasing trend of δ13C in fish collected in the Sea of Japan and East China Sea (Ohshimo et al., 2021) was supported by the decreasing trends of δ13C in their prey and is possibly explained by the Suess effect. Therefore, the changing carbon environments are recorded in the marine ecosystem in the Sea of Japan and are detectable by the stable isotope approach even up to 15 years in the past.

The interannual linear trend of the residual δ15N was not significant. In the South China Sea and other western North Pacific marginal seas, anthropogenic changes in nitrogen dynamics are detectable using the stable isotope approach (Ren et al., 2017). These results suggest that nitrogen dynamics in the Sea of Japan do not change rapidly compared to those in the South China Sea, although different organisms are measured. Our results also differed from the trend of small pelagic fish in the Sea of Japan and East China Sea reported by Ohshimo et al. (2021). A decreasing trend of δ15N ( 0.05 ‰ yr−1) was detected in small pelagic fish in the Sea of Japan and East China Sea based on a time series analysis. A significant decreasing trend of δ15N was not detected in the GLM approach, but low δ15N values in small pelagic fish are recorded in 2016 and 2017 (after 2005) (Ohshimo et al., 2021). The residual δ15N values in this study were not very low in these 2 years, suggesting that the δ15N values of secondary producers and small pelagic fish were not coupled in our study area. Additionally, some studies have reported eventual phenomena in the Sea of Japan. For example, considerable amounts of nitrate were discharged into the Sea of Japan by the Changjiang River during the summer of 2013 (Kodama et al., 2017b), and the surface salinity was lowered by inflow from this river in the summer or autumn of 2010, 2012, and 2015 (Kosugi et al., 2021). However, we could not find any relationship between the occasional change in the stable isotope ratios of C. sinicus and the reported eventual phenomena in the Sea of Japan. These results suggest that nitrogen dynamics supporting biological production in this area are very complex, and long-term monitoring is required to identify the anthropogenic impact on nitrogen dynamics in the Sea of Japan.

5 Conclusions

We used data from a 15-year study of δ13Cbulk and δ15Nbulk values of the calanoid copepod C. sinicus to indirectly examine the interannual variations of isotope ratios of stable carbon and nitrogen. The δ15Nbulk values confirmed that C. sinicus is a “secondary” producer in the area. The δ13Cbulk values of C. sinicus reflected both the δ13C of phytoplankton and the physiological processes of C. sinicus, and the variations in δ13Cbulk were largely explained by environmental parameters, geographical position, and elemental ratios of C. sinicus using the GLM approach. The residuals between the observed stable isotope ratios and GLM-predicted stable isotope ratios (residual δ13C and residual δ15N) indicated that δ13C showed a significant linear decreasing trend, but this was not observed for residual δ15N. The decreasing trend of residual δ13C was 0.0365 ‰ yr−1± 0.0124 ‰ yr−1, which does not differ considerably from that observed in the muscle of small pelagic fish in the East China Sea and the Sea of Japan owing to the Suess effect. Therefore, anthropogenic effects on carbon dynamics were recorded and detectable in the Sea of Japan even during the 15-year monitoring period. However, long-term monitoring is required to identify the anthropogenic impact on nitrogen dynamics in the Sea of Japan. The δ15N results also indicated that nitrogen dynamics supporting the biological production are complex in the Sea of Japan, and more investigations are required to understand the overview of the nitrogen dynamics in this area, such as the seasonal variations in the δ15N of POM.

Code and data availability

Our code and data are in the Mendeley Data repository (Nakamura et al., 2021;


The supplement related to this article is available online at:

Author contributions

KN, NI, HM, and TK designed the experiments. KN, AN, SY, YK, MN, NI, HM, and TK carried them out. KN, AN, and TK prepared the paper with contributions from all co-authors.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Onboard observations were conducted through R/Vs Mizuho-Maru and Yoko-Maru of the Japan Fisheries Research and Education Agency as well as Dai-Roku Kaiyo-Maru of Kaiyo Engineering Co., Ltd. We appreciate captains, crews, researchers, and staff supporting us and participating in sampling during the cruises. We also thank Keiko Yamada and Kyoko Matsuda for support of onshore experiments.

Financial support

This research has been supported by the Japan Society for the Promotion of Science (grant nos. 16K07831 and 19K06198), the Fisheries Agency of Japan, and the Japan Fisheries Research and Education Agency.

Review statement

This paper was edited by Yuelu Jiang and reviewed by four anonymous referees.


Aita, M. N., Tadokoro, K., Ogawa, N. O., Hyodo, F., Ishii, R., Smith, S. L., Saino, T., Kishi, M. J., Saitoh, S.-I., and Wada, E.: Linear relationship between carbon and nitrogen isotope ratios along simple food chains in marine environments, J. Plankton Res., 33, 1629–1642,, 2011. 

Antonio, E. S., Kasai, A., Ueno, M., Ishihi, Y., Yokoyama, H., and Yamashita, Y.: Spatial-temporal feeding dynamics of benthic communities in an estuary-marine gradient, Estuar. Coast. Shelf. Sci., 112, 86–97,, 2012. 

Chen, C.-T. A., Lui, H.-K., Hsieh, C.-H., Yanagi, T., Kosugi, N., Ishii, M., and Gong, G.-C.: Deep oceans may acidify faster than anticipated due to global warming, Nat. Clim. Change, 7, 890–894,, 2017. 

Doney, S. C.: The growing human footprint on coastal and open-ocean biogeochemistry, Science, 328, 1512–1516,, 2010. 

Duce, R. A., LaRoche, J., Altieri, K., Arrigo, K. R., Baker, A. R., Capone, D. G., Cornell, S., Dentener, F., Galloway, J., Ganeshram, R. S., Geider, R. J., Jickells, T., Kuypers, M. M., Langlois, R., Liss, P. S., Liu, S. M., Middelburg, J. J., Moore, C. M., Nickovic, S., Oschlies, A., Pedersen, T., Prospero, J., Schlitzer, R., Seitzinger, S., Sorensen, L. L., Uematsu, M., Ulloa, O., Voss, M., Ward, B., and Zamora, L.: Impacts of atmospheric anthropogenic nitrogen on the open ocean, Science, 320, 893–897,, 2008. 

Furuichi, S., Yasuda, T., Kurota, H., Yoda, M., Suzuki, K., Takahashi, M., and Fukuwaka, M.: Disentangling the effects of climate and density-dependent factors on spatiotemporal dynamics of Japanese sardine spawning, Mar. Ecol. Prog. Ser., 633, 157–168,, 2020. 

Gatten, R. R., Sargent, J. R., Forsberg, T. E. V., O'Hara, S. C. M., and Corner, E. D. S.: On the nutrition and metabolism of zooplankton XIV. Utilization of lipid by Calanus helgolandicus during maturation and reproduction, J. Mar. Biol. Assoc. UK, 60, 391–399,, 1980. 

Gruber, N., Keeling, C. D., Bacastow, R. B., Guenther, P. R., Lueker, T. J., Wahlen, M., Meijer, H. A. J., Mook, W. G., and Stocker, T. F.: Spatiotemporal patterns of carbon-13 in the global surface oceans and the oceanic suess effect, Global Biogeochem. Cy., 13, 307–335,, 1999. 

Halpern, B. S., Walbridge, S., Selkoe, K. A., Kappel, C. V., Micheli, F., D'Agrosa, C., Bruno, J. F., Casey, K. S., Ebert, C., Fox, H. E., Fujita, R., Heinemann, D., Lenihan, H. S., Madin, E. M., Perry, M. T., Selig, E. R., Spalding, M., Steneck, R., and Watson, R.: A global map of human impact on marine ecosystems, Science, 319, 948–952,, 2008. 

Hirai, J., Hamamoto, Y., Honda, D., and Hidaka, K.: Possible aplanochytrid (Labyrinthulea) prey detected using 18S metagenetic diet analysis in the key copepod species Calanus sinicus in the coastal waters of the subtropical western North Pacific, Plank. Benth. Res., 13, 75–82,, 2018. 

Hirakawa, K. and Goto, T.: Diet of larval sardine, Sardinops melanostictus in Toyama Bay, southen Japan Sea, Bull. Japan Sea Reg. Fish. Res. Lab., 46, 65–75, 1996. 

Hirakawa, K., Goto, T., and Hirai, M.: Diet composition and prey size of larval anchovy, Eugraulis japonicus, in Toyama Bay, Southern Japan Sea, Bull. Japan Sea Reg. Fish. Res. Lab., 47, 67–78, 1997. 

Holm-Hansen, O., Lorenzen, C. J., Holmes, R. W., and Strickland, J. D.: Fluorometric determination of chlorophyll, J. Conseil, 30, 3–15, 1965. 

Huang, C., Uye, S., and Onbé, T.: Geographic distribution, seasonal life cycle, biomass and production of a planktonic copepod Calanus sinicus in the Inland Sea of Japan and its neighboring Pacific Ocean, J. Plankton Res., 15, 1229–1246,, 1993. 

Ishizu, M., Miyazawa, Y., Tsunoda, T., and Ono, T.: Long-term trends in pH in Japanese coastal seawater, Biogeosciences, 16, 4747–4763,, 2019. 

Jager, T., Salaberria, I., and Hansen, B. H.: Capturing the life history of the marine copepod Calanus sinicus into a generic bioenergetics framework, Ecol. Model., 299, 114–120,, 2015. 

Kitayama, K., Seto, S., Sato, M., and Hara, H.: Increases of wet deposition at remote sites in Japan from 1991 to 2009, J. Atmos. Chem., 69, 33–46,, 2012. 

Kodama, T., Morimoto, H., Igeta, Y., and Ichikawa, T.: Macroscale-wide nutrient inversions in the subsurface layer of the Japan Sea during summer, J. Geophys. Res., 120, 7476–7492,, 2015. 

Kodama, T., Igeta, Y., Kuga, M., and Abe, S.: Long-term decrease in phosphate concentrations in the surface layer of the southern Japan Sea, J. Geophys. Res., 121, 7845–7856,, 2016. 

Kodama, T., Hirai, J., Tamura, S., Takahashi, T., Tanaka, Y., Ishihara, T., Tawa, A., Morimoto, H., and Ohshimo, S.: Diet composition and feeding habits of larval Pacific bluefin tuna Thunnus orientalis in the Sea of Japan: integrated morphological and metagenetic analysis, Mar. Ecol. Prog. Ser., 583, 211–226,, 2017a. 

Kodama, T., Morimoto, A., Takikawa, T., Ito, M., Igeta, Y., Abe, S., Fukudome, K., Honda, N., and Katoh, O.: Presence of high nitrate to phosphate ratio subsurface water in the Tsushima Strait during summer, J. Oceanogr., 73, 759–769,, 2017b. 

Kodama, T., Wagawa, T., Iguchi, N., Takada, Y., Takahashi, T., Fukudome, K. I., Morimoto, H., and Goto, T.: Spatial variations in zooplankton community structure along the Japanese coastline in the Japan Sea: influence of the coastal current, Ocean Sci., 14, 355–369,, 2018a. 

Kodama, T., Wagawa, T., Ohshimo, S., Morimoto, H., Iguchi, N., Fukudome, K. I., Goto, T., Takahashi, M., and Yasuda, T.: Improvement in recruitment of Japanese sardine with delays of the spring phytoplankton bloom in the Sea of Japan, Fish. Oceanogr., 27, 289–301,, 2018b. 

Kodama, T., Nishimoto, A., Horii, S., Ito, D., Yamaguchi, T., Hidaka, K., Setou, T., and Ono, T.: Spatial and seasonal variations of stable isotope ratios of particulate organic carbon and nitrogen in the surface water of the Kuroshio, J. Geophys. Res., 126, e2021JC017175,, 2021. 

Kogure, Y.: Stable carbon and nitrogen isotope analysis of the sublittoral benthic food web structure of an exposed sandy beach, Bull. Biogeogr. Soc. Jpn., 59, 15–25, 2004. 

Kosugi, N., Hirose, N., Toyoda, T., and Ishii, M.: Rapid freshening of Japan Sea Intermediate Water in the 2010s, J. Oceanogr., 77, 269–281,, 2021. 

Logan, J. M., Jardine, T. D., Miller, T. J., Bunn, S. E., Cunjak, R. A., and Lutcavage, M. E.: Lipid corrections in carbon and nitrogen stable isotope analyses: comparison of chemical extraction and modelling methods, J. Anim. Ecol., 77, 838–846,, 2008. 

Lorrain, A., Pethybridge, H., Cassar, N., Receveur, A., Allain, V., Bodin, N., Bopp, L., Choy, C. A., Duffy, L., Fry, B., Goni, N., Graham, B. S., Hobday, A. J., Logan, J. M., Menard, F., Menkes, C. E., Olson, R. J., Pagendam, D. E., Point, D., Revill, A. T., Somes, C. J., and Young, J. W.: Trends in tuna carbon isotopes suggest global changes in pelagic phytoplankton communities, Global Change Biol., 26, 458–470,, 2020. 

Lüdecke, D.: ggeffects: Tidy Data Frames of Marginal Effects from Regression Models, J. Open Source Softw., 3, 772,, 2018. 

Mayor, D. J., Cook, K., Thornton, B., Walsham, P., Witte, U. F. M., Zuur, A. F., and Anderson, T. R.: Absorption efficiencies and basal turnover of C, N and fatty acids in a marine Calanoid copepod, Funct. Ecol., 25, 509–518,, 2011. 

Nakamura, K.-i., Nishimoto, A., Yasui-Tamura, S., Kogure, Y., Nakae, M., Iguchi, N., Morimoto, H., and Kodama, T.: Data for: “Stable isotope ratio of Calanus sinicus in the coastal Japan Sea” (2) [data set],, 2021. 

Nishida, K., Yasu, A., Nanjo, N., Takahashi, M., Kitajima, S., and Ishimura, T.: Microscale stable carbon and oxygen isotope measurement of individual otoliths of larvae and juveniles of Japanese anchovy and sardine, Estuar. Coast. Shelf. Sci., 245, 106946,, 2020. 

Ohman, M. D., Rau, G. H., and Hull, P. M.: Multi-decadal variations in stable N isotopes of California Current zooplankton, Deep Sea Res. I, 60, 46–55,, 2012. 

Ohshimo, S., Tawa, A., Ota, T., Nishimoto, S., Ishihara, T., Watai, M., Satoh, K., Tanabe, T., and Abe, O.: Horizontal distribution and habitat of Pacific bluefin tuna, Thunnus orientalis, larvae in the waters around Japan, Bull. Mar. Sci., 93, 769–787,, 2017. 

Ohshimo, S., Madigan, D. J., Kodama, T., Tanaka, H., Komoto, K., Suyama, S., Ono, T., and Yamakawa, T.: Isoscapes reveal patterns of δ13C and δ15N of pelagic forage fish and squid in the Northwest Pacific Ocean, Prog. Oceanogr., 175, 124–138,, 2019. 

Ohshimo, S., Kodama, T., Yasuda, T., Kitajima, S., Tsuji, T., Kidokoro, H., and Tanaka, H.: Potential fluctuation of δ13C and δ15N values of small pelagic forage fish in the Sea of Japan and East China Sea, Mar. Freshwater Res., 72, 1811,, 2021. 

Ono, T.: Long-term trends of oxygen concentration in the waters in bank and shelves of the Southern Japan Sea, J. Oceanogr., 77, 659-684,, 2021. 

Pu, X.-M., Sun, S., Yang, B., Zhang, G.-T., and Zhang, F.: Life history strategies of Calanus sinicus in the southern Yellow Sea in summer, J. Plankton Res., 26, 1059–1068,, 2004a. 

Pu, X.-M., Sun, S., Yang, B., Ji, P., Zhang, Y.-S., and Zhang, F.: The combined effects of temperature and food supply on Calanus sinicus in the southern Yellow Sea in summer, J. Plankton Res., 26, 1049–1057,, 2004b. 

R Core Team: R: A language and environment for statistical computing, [code], 2022.  

Ren, H., Chen, Y. C., Wang, X. T., Wong, G. T. F., Cohen, A. L., DeCarlo, T. M., Weigand, M. A., Mii, H. S., and Sigman, D. M.: 21st-century rise in anthropogenic nitrogen deposition on a remote coral reef, Science, 356, 749–752,, 2017. 

Smyntek, P. M., Teece, M. A., Schulz, K. L., and Thackeray, S. J.: A standard protocol for stable isotope analysis of zooplankton in aquatic food web research using mass balance correction models, Limnol. Oceanogr., 52, 2135–2146,, 2007. 

Tanaka, H., Takasuka, A., Aoki, I., and Ohshimo, S.: Geographical variations in the trophic ecology of Japanese anchovy, Engraulis japonicus, inferred from carbon and nitrogen stable isotope ratios, Mar. Biol., 154, 557–568,, 2008. 

Uye, S.: Why does Calanus sinicus prosper in the shelf ecosystem of the Northwest Pacific Ocean?, ICES J. Mar. Sci., 57, 1850–1855,, 2000. 

Uye, S. and Yamamoto, F.: In situ feeding of the planktonic copepod Calanus sinicus in the Inland Sea examined by the gut fluorescence method, Bull. Plankton Soc. Japan, 42, 123–139, 1995. 

Uye, S.-i.: Temperature-dependent development and growth of Calanus sinicus (Copepoda: Calanoida) in the laboratory, Hydrobiologia, 167, 285–293,, 1988. 

Yi, X., Huang, Y., Zhuang, Y., Chen, H., Yang, F., Wang, W., Xu, D., Liu, G., and Zhang, H.: In situ diet of the copepod Calanus sinicus in coastal waters of the South Yellow Sea and the Bohai Sea, Acta Oceanol. Sin., 36, 68–79,, 2017. 

Zhang, G.-T., Sun, S., and Zhang, F.: Seasonal variation of reproduction rates and body size of Calanus sinicus in the Southern Yellow Sea, China, J. Plankton Res., 27, 135–143,, 2005. 

Zhou, K., Sun, S., Wang, M., Wang, S., and Li, C.: Differences in the physiological processes of Calanus sinicus inside and outside the Yellow Sea Cold Water Mass, J. Plankton Res., 38, 551–563,, 2016. 

Short summary
The Sea of Japan, surrounding Russia, the Korean Peninsula, and the Japanese Archipelago, is one of the most rapidly changing seas in the world. We measured carbon and nitrogen isotope ratios of zooplankton. We determined that the carbon stable isotope ratio has been decreasing over 15 years, and this trend was comparable to or slightly more rapid than the Suess effect, which is a signal of anthropogenic disturbance. Therefore, carbon dynamics are changing in the shallow coastal waters.