All statistical analyses were conducted using R (R Core Team, 2023). Accompanied by the analysis of variance (ANOVA), a pairwise test with the Tukey–Kramer adjacent was applied to the least-squared mean (lsmean) values. The “ggeffect” package (Lüdecke, 2018) was used to calculate the lsmean values and SEs. The “car” package (Fox and Weisberg, 2018) of ANOVA was used to conduct type II ANOVA for unbalanced data comparison. Neither δ13CPOM nor δ15NPOM exhibited a normal distribution according to the Kolmogorov–Smirnov test (p < 0.001); both showed multimodal distributions. Therefore, we applied a two-dimensional Gaussian mixed model (GMM) to classify δ13CPOM and δ15NPOM using the “mclust” package (Scrucca et al., 2016). The samples were not divided based on the sampling depth and seasons. The number of classes was determined using the Bayesian information criterion (BIC).