1 Introduction

Cephalopods (e.g., squid, octopod, and cuttlefish) are a major component of global marine ecosystems (Chesnais et al., 2019). A wide variety of species representing important components of the local ecosystem (Clarke, 1996; Xavier et al., 2014) inhabit coastal and shelf waters, which are characterized by high primary and secondary productivities (Longhurst et al., 1995). Cephalopods in coastal and oceanic areas have recently increased in abundance in response to changing ocean environments (Doubleday et al., 2016; Pang et al., 2018). Ce phalopods are characterized by rapid growth, a short lifespan, high turnover rates, and strong phenotypic plasticity (Pecl and Jackson, 2008; Arkhipkin et al., 2021). Their flexible life strategies render them capable of quick responses to changing environments. Researchers believe that cephalopods may eventually overtake fish species in coastal ecosystems, which have become degraded as a consequence of fishing pressure and climate change (Pierce et al., 2008; Doubleday et al., 2016).

China's seas, which are among the largest marine ecosystems in the world, have demonstrated rapid warming with climate change (Pang et al., 2018). For example, sea surface temperatures (SSTs) in the East China Sea have risen by 1.22℃ from 1982 to 2006 (Belkin, 2009); in the Yellow Sea (YS), the SST has increased by 1.7℃ since 1980 (Lin et al., 2005; Pei et al., 2017). China's seas are overexploited, and commercial fish landings have decreased significantly; interestingly, however, cephalopod landings have increased (Pang et al., 2018). The YS is a semi-enclosed marginal sea located in the Northwest Pacific Ocean, surrounded by the Korean Peninsula and mainland China. The area of this sea is approximately 380000 km², and its average depth is approximately 44 m (Park et al., 2015). The YS is an intensively exploited shallow sea ecosystem (Wang et al., 2016) and an important global resource for coastal and offshore fisheries (Tang, 2003). Climate change, overfishing, and exploitation have affected the self-regulating mechanism of the YS ecosystem (Tang, 2003). Indeed, 90% of its fishery resources are now in danger of exhaustion change (Jin, 2017), which could alter its community structure (Lin et al., 2005; Wang et al., 2016). Despite this situation, the cephalopod catch in this region has increased since 2003, stabilizing at approximately 10 × 10⁴ t from 2003 to 2013 (Wang et al., 2016). Chlorophyll a (Chl-a) concentrations have also trended upward since 1998 (Tan et al., 2011; Tian et al., 2020).

The cephalopod fauna of the YS mainly comprises 14 warm-temperate species belonging to six families. The dominant taxa are Loliolus and Octopus (mainly L. beka and A. fangsiao), followed by Todarodes pacificus (Du et al., 2016). L. beka and A. fangsiao are the primary catch species of cephalopod fisheries, and they are commercially and ecologically important throughout the whole China's waters (Dong, 1978; Du et al., 2016; Jin et al., 2017). For example, the A. fangsiao catch in early spring of 1970s was 85% of the annual cephalopod catch in China (Dong, 2014). Both species live in shallow coastal waters and have life cycles of less than one year. The loliginid L. beka performs diel migrations in surface waters at night, while postparalarva A. fangsiao is basically benthic (Dong, 1988; Jin et al., 2020; Pang et al., 2020). Both species are highly sensitive to environmental changes, which renders them of potential value in monitoring exercises to detect environmental change (Pecl and Jackson, 2008; Chesnais et al., 2019). As one species is a squid and the other is an octopus, they also likely reflect different responses to environmental changes. Research on these taxa has focused on molecular genetic analysis (Jiang et al., 2015; Um et al., 2017). Unfortunately, the lack of basic biological information on these species limits our acquisition of a comprehensive understanding of their population dynamics.

Condition factors are important for evaluating the health of aquatic species, populations, and communities (LeCren, 1951; Ferreri, 2014). They can be affected by an individual's size, sex, sexual cycle, food supply, and habitat (Mohapatra et al., 2010; Zhang et al., 2017). Without specific age data, condition factors have indicated biochemical, ecological, and physiological processes in studies on the phenotypic responses of fish (Maina et al., 2019), crab (Zhang et al., 2017), shrimp (Ramírez and Paramo, 2020), and cephalopods (Ferreri, 2014; Siddique et al., 2014) to environmental change. Growth integrates the effects of physical and biological processes and is an ideal parameter for measuring the response of organisms to environmental change and fishing pressure (Morrongiello et al., 2012).

The scarcity of information on the population dynamics of cephalopods in the YS hinders our understanding of how cephalopods have already responded and will respond to changes in the marine environment. In this study, we investigate the appropriateness of L. beka and A. fangsiao as indicator species for detecting environmental change by comparing inter-annual variations in their body weight (BW). We use a mixed-effects model to analyze the weight growth of these two species responding to environmental factors with the goals of 1) reproducing inter-annual variations in BW for both species, and 2) analyzing environmental factors on the weight growth of each species.

2 Materials and Methods 2.1 Sampling Survey

Squid and octopus samples were collected from bottom trawl surveys in the YS during the spring (May) and autumn (September) of 2011 and 2013 – 2018 (Fig.1). The survey vessel was a 220 kW single tugboat that was carried out during the day with a towing speed of 2 – 3 knots per hour at each station. We determined the water depth at each station through measurements of conductivity, temperature, and depth profiler and recorded each tow. The survey water depth ranged from 10 m to 75 m, and the average water depth was 35.2 m. Species were identified following Dong (1988) and Okutani (2015). The dorsal mantle length (ML, mm) and BW (0.1 g) of 4192 L. beka and 1100 A. fangsiao individuals (representing 66.21% and 19.32% of the total cephalopod catch, respectively; Table 1) were measured.

2.2 Marine Environmental Data

SST and Chl-a are oceanographic variables that are considered most suitable for determining cephalopod habitat (Arkhipkin et al., 2015). The Southern Oscillation Index (SOI) was used as an important climate index in the Pacific (Power and Kociuba, 2011). SOI is significantly negatively correlated with the average SST in the YS and affects the intensity of the YS warm current, thereby impacting the YS ecosystem (Yin et al., 2016). We selected SST, Chl-a, and SOI as environmental variables to study the effects of their change on the BW of L. beka and A. fangsiao over time. Data for each species from 2010 to 2018 were sourced as follows. (a) Monthly average SST and Chl-a data (119˚ – 124˚E, 34˚ – 37.5˚N) were derived from the MODIS (Moderate-Resolution Imaging Spectroradiometer) Aqua products of the National Oceanic and Atmospheric Administration (NOAA) (https://coastwatch.pfeg.noaa.gov/). (b) Monthly average SOI data were derived from the NOAA (https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/soi.long.data). Here, the annual average SST and Chl-a are represented by the annual average of their corresponding monthly averages (± SE). The spatial resolution of all marine environmental data is 0.1˚.

2.3 Statistical Analysis and Model Building

BW is used to evaluate the condition factors of L. beka and A. fangsiao. Regression and correlation analyses were performed between ML and BW, and a one-way analysis of variance was used to identify differences in weight growth among years. A least-squares method was used to test length-weight relationships (LWR) for both species with following equation (LeCren, 1951):

$ BW = aM{L^b}, $

(1)

where a represents the individual's nutritional status, and b represents its growth pattern. Because the variance of BW increases as ML increases, the above equation was logtransformed to:

$ \ln (BW) = \ln (a) + b \times \ln (ML). $

(2)

We fitted a generalized linear model (GLM), a simple linear model suitable for individuals at different years and water depths (for L. beka and A. fangsiao, respectively), and two linear mixed-effect models (LMEM) considering the year and/or depth as random factors to explain the relationship between BW and ML (Ma et al., 2017). All models for ML included the maximum intrinsic fixed effect structure. Akaike's information criterion for small sample size (AICc) was used to evaluate the relative support of each candidate set in the model, which was further corrected for the small sample size, to determine the optimal growth model (Beier et al., 2001) (Table 2).

Sliding window analysis was performed on the basis of the growth model to determine the optimal time window for weight growth to occur in response to environmental variables (Van de Pol et al., 2016). This type of analysis can test the influence of different environmental variables on the response variables in different time windows and compare different time windows to determine which environmental variables can best predict biological responses (Smoliński, 2019b). Absolute time windows and linear and quadratic relationships were assumed in the analysis. In addition, 1000 randomization tests were performed to eliminate the possible misclassification of environmental signals due to change (i.e., false positives) (Smoliński, 2019a). The possibility of obtaining a powerful model supported by chance was quantified by calculating the P_∆AICc value from 1000 randomized trials. The sliding window analysis independently considers the external effects of SST, Chl-a, and SOI as fixed effects within a certain time window. The AICc and significance tests were used to determine the best model, and restricted maximum likelihood estimation was used for refitting to obtain unbiased variable estimates (Zuur et al., 2010). The predicted effect P-value shows the degree of correlation between environmental variables and BW, and the t value represents the sequential t-test analysis (STARS) results of the obtained time series to identify significant step changes in BW. All analyses were performed using the R Studio 3.5; 'lme4' was used for building mixed models; 'MuMIn' was used to calculate conditional R²-values; 'AICcmodavg' was used to calculate AICc; and 'climwin' was used to obtain the sliding window to predict which environmental variables affect weight growth the most.

3 Results 3.1 Mantle Length and Body Weight Distribution

Inter-annual changes in the frequency distributions of ML and BW in spring and autumn are apparent for L. beka and A. fangsiao (Fig.2). The ML distributions of both species in autumn are generally unimodal (mainly between 30 and 40 mm for L. beka, and between 40 and 50 mm for A. fangsiao), whereas ML distributions during spring are more scattered. The BW distribution of L. beka is concentrated between 0 and 5 g in spring and autumn, whereas that of A. fangsiao is relatively dispersed in the same seasons. The range and mean of the BW of A. fangsiao in spring are higher than those in autumn.

3.2 Length-Weight Relationship

The LWR relationships (as Eq. (1)) of L. beka and A. fangsiao differ seasonally (Fig.3).

For L. beka,

$ {W_{{\text{spring}}}} = 0.0015M{L^{2.1243}} (n = 1496, {R^2} = 0.7632), $

$ {W_{{\text{autumn}}}} = 0.0036M{L^{1.8906}} (n = 2696, {R^2} = 0.6258). $

For A. fangsiao,

$ {W_{{\text{spring}}}} = 0.0103M{L^{2.1446}} (n = 264, {R^2} = 0.7918), $

$ {W_{{\text{autumn}}}} = 0.0143M{L^{2.0348}} (n = 836, {R^2} = 0.6968). $

Whereas the allometric coefficient b for each species is the largest in spring, a is the largest in autumn. In general, BW and ML values in autumn are lower than those in spring, which indicates that spring environments are more conducive to growth.

3.3 Inter-Annual Changes in Weight Growth

AICc values indicate that LMEM with year and depth and random effects on the intercept is the best model (Table 2) among those tested, thus reflecting the impact of environmental differences in the habitat depth of these two species on their BW. BW models for L. beka and A. fangsiao respectively explain 73.80% and 75.26% of the variance noted for each species, and the 7-year weight growths of L. beka and A. fangsiao are reconstructed based on the growth model. The BW of each species varies greatly among different years, and the trends observed vary between species (Fig.4). The BW of L. beka manifests an overall upward trend, especially after 2014, and peaks in 2017. The BW of A. fangsiao fluctuates remarkably and is the lowest in 2013, and the highest in 2016. Obvious seasonal differences in intraspecific BW patterns are apparent (Fig.5). The maximum peak for A. fangsiao occurs in spring and autumn of 2016, while that for L. beka occurs in spring of 2017. The growth trends gradually increase from autumn of 2011. Both annual (Fig.4) and seasonal (Fig.5) BW changes show obvious annual variations. The 'year' factor plays a key role in understanding the long-term series BW of the two species, while the 'season' factor exerts a local effect on the current year.

3.4 Inter-Annual Variations in Environmental Factors

Inter-annual variations in SST, Chl-a, and SOI are shown in Fig.6. SST increases from 13.8℃ in 2011 to 16.6℃ in 2016, and then to 16.8℃ in 2017. Chl-a trends downward from 2014 (2.8 mg m⁻³) and is the lowest in 2018 (2.0 mg m⁻³). SOI decreases from 2011 to 2015, increases from 2015, and peaks in 2011 (1.4). The lowest SOIs are noted in 2015 (−1.2) and 2016 (−0.5).

3.5 Environmental Effects on Body Weight

Sliding window analysis was used to identify the optimal time window for each environmental variable (Table 3; Figs.7, 8). The inclusion of these variables improves the model fitness and is supported by the AICc comparisons (Figs.7, 8). The ΔAICc value of the best model fitted to the observed data was compared with that of the null model to obtain the optimal opening and closing time window periods, and the environmental variables all reveal strong signals (red areas in Figs.7, 8). The BW response of both species is significantly related to SST and SOI, and the P_∆AICc obtained from the randomization test is < 0.001. These results suggest that SST and SOI signals are very unlikely to be false positives. In a different response window, BW increases with increasing SST and decreases with increasing SOI (Figs.9, 10). L. beka also responds significantly to Chl-a, and the P_∆AICc determined from the randomization test is < 0.001. This value indicates that the signal is unlikely to be a false positive (Fig.9).

4 Discussion

Because of rapid warming and overexploitation of China's seas, the YS ecosystem has degraded, and high-trophic level species have decreased (Tang, 2003). However, the abundance of cephalopods has increased, most likely because of their significant phenotypic plasticity and reduction in the number of predators (Wang et al., 2016; Arkhipkin et al., 2021). L. beka and A. fangsiao occupy relatively high trophic levels and feed mainly on benthic crustaceans and small fish larvae (Huang, 2004; FAO, 2010). These two species have come to occupy dominant positions in the YS ecosystem (Du et al., 2016).

4.1 Growth Characteristics

The developments of cephalopods differ as a result of differences in their life history, gonad development, food supplies, and habitats (Semmens et al., 2004; Pecl and Jackson, 2008; Pang et al., 2020). We found that the BW of L. beka and A. fangsiao varies remarkably over time, and the trends for each species differ (Fig.4) according to their unique life histories (Dong, 1988). Obvious seasonal differences in the growth of the same species can be observed (Fig.5), which is closely related to their changing life history characteristics at different developmental stages and seasonality in their habitat (Semmens et al., 2004; Pecl and Jackson, 2008). The number of individuals of each species caught in spring is lower than that in autumn (Table 1), which may be related to fishery seasonality (Jin et al., 2020). The spawning population of L. beka is mainly targeted in spring, and lower catches are often obtained in this season; individuals weighing < 5 g account for 69.45% of the total catch in this season (Fig.2). A. fangsiao is mainly obtained in early spring when it begins to spawn; the catch during this season is usually smaller than that in other seasons.

We could only investigate seasonal differences in LWR because many individuals in the samples were unsexed (Fig.3). LWR differs between species and seasonally within a species, as mainly reflected by parameters a and b (Fig.3). Because the LWR b values of the two species are < 3, their weight growth follows an allometric growth pattern. Moreover, the b values of both species are greater in spring than those in autumn, indicating greater spring BWs. This finding is supported by the data in Table 1. The numbers and weights of individuals of ML > 60 mm in spring are greater than those in autumn (Fig.3), which could be related to overfishing after the closure of the summer fishing closure season. The LWR R² values for each species are relatively low because growth is susceptible to changes in environmental conditions (Semmens et al., 2004). This result is consistent with other studies on octopus species (Silva et al., 2002; Otero et al., 2007; Pang et al., 2020). Previous studies showed that A. fangsiao mainly spawns in spring (Dong, 2014), but the spawning season of the species may actually range from early spring to late autumn (Pang et al., 2020). A few of studies have reported the growth of L. beka. Large variations in ML and BW are also noted between individuals in the same season (Table 1, Fig.3), which may be caused by different developmental stages and age structures, as well as other factors. Further researches combining age determination and growth pattern analysis are recommended on both species.

4.2 Effects of Environmental Factors on Growth

Temperature affects the aggregation, reproduction, migration, growth, abundance, and population dynamics of cephalopods (Arkhipkin et al., 2015; Angeles-Gonzalez et al., 2017). L. beka and A. fangsiao demonstrate a positive linear growth-SST relationship (Figs.9, 10) and similar BW patterns from 2014 – 2016 and during spring (Figs. 4, 5). Similar positive linear respond patterns have been observed in L. forbesi, L. pealeii, Illex coindetii, Uroteuthis duvaucelii (Forsythe, 2004; Keller et al., 2017; Sasikumar et al., 2018). Both species inhabit shallow waters, where temperature can synchronize life history features in seasonal habitats (Quetglas et al., 2011). In addition, the patterns of change in the inter-annual BW of L. beka roughly mirror those of SST (Figs.4, 6). BW is the highest in 2017 and the lowest in 2011 (Fig.4), mainly because warmer waters support higher metabolic rates and higher food intake (Semmens et al., 2004). Compared with the results of other studies in the past 20 years, cephalopods in the YS appear to have migrated further north because of warmer temperatures (Jin et al., 2020). Moreover, key window periods for L. beka and A. fangsiao growths occur in June and May – August (Table 3) when temperatures are 14.9 – 22.6℃ and 17.5 – 23.4℃, respectively. These temperatures differ from those reported in previous studies (Jin et al., 2020), possibly because of differences in seasonal habitats. Further research is needed to examine the impacts of climate change on seasonal species to predict changes in the YS ecosystem.

Chl-a strongly influences prey density, habitat suitability, catch per unit effort (CPUE), and cephalopod abundance (Puerta et al., 2015; Keller et al., 2017; Yu et al., 2018). We observed a strong negative correlation between L. beka weight and Chl-a (Fig.9), possibly because L. beka performs diel migrations near the seabed during the day and near the surface at night. By contrast, A. fangsiao is a benthic species, and its response to Chl-a is delayed by several months (Dong, 1988; Puerta et al., 2015; Jin et al., 2020). Inter-annual variations in L. beka BW trend in the opposite direction to those of Chl-a in our study. While Chl-a trends upward from 2011 to 2014 and downward from 2014 to 2018 (Fig.6), the BW of L. beka is poor from 2011 – 2014 and improves after 2014 (Fig.4). Chl-a in the corresponding critical window (mainly during winter and spring) is in the high-value during the year (Tian et al., 2020), but the BW of L. beka is poor, which may be related to interspecies competition (Puerta et al., 2015; Keller et al., 2017). Because small pelagic fishes effectively compete for food with early and juvenile stages of L. beka, higher Chl-a may promote the growth of the former, thereby reducing the prey biomass for L. beka, which leads to its reduced weight (Keller et al., 2017). However, the negative correlation between Chl-a and growth may also be related to the growth rate. Higher primary productivity can accelerate development, thereby hastening sexual maturity at a smaller body size (Semmens et al., 2004). The Chl-a range (2.4 – 3.4 mg m⁻³) in this critical window is within the appropriate range (2.3 – 3.76 mg m⁻³) reported in previous studies, but the relationship between Chl-a and BW is contradictory (Jin et al., 2020). This observation may be related to habitat seasonality and reflects the impact of environmental changes on species adaptability.

Climate change affects the growth, reproduction, and behavior of cephalopods, especially during the El Niño-Southern Oscillation (ENSO) (Yu et al., 2018). ENSO events, especially the SST anomaly and intensity of warm currents in winter, affect SST by changing the surface heat flux within the YS; the strength of the main flow of the warm current has a lag relationship with the SOI (Wang et al., 2020). We found that the BW and SOI of these two species have a significant negative linear relationship, with the sliding window in the El Niño period exerting a delayed effect on growth, consistent with previous studies (Chen et al., 2007; Yu et al., 2018; Wang et al., 2020). The consistent response of the BW of the studied species to the SOI confirms that ENSO can synchronize aspects of the life histories of a wide range of taxa over large geographic scales (Black et al., 2018; Tanner et al., 2020). El Niño and La Niña events affect the biological and physical environments, such as the hatching and feeding conditions of spawning grounds, thereby potentially affecting the number and growth of individuals likely to recruit into subsequent populations (Chen et al., 2007). Because El Niño was strong during 2015 – 2016, temperatures in the YS increased, and the BW of L. beka and A. fangsiao increased. By contrast, 2011 was a La Niña year, and the lower temperatures in the YS were not conducive to the BW of both species (Figs.4, 6). Studies have shown that the spring and autumn CPUE of A. fangsiao in Haizhou Bay may be significantly negatively correlated with temperature during El Niño periods (Pang et al., 2020). ENSO events also affect other fish and zooplankton communities in the YS (Liu et al., 2020; Shi et al., 2020). For example, increased SST caused the wintering grounds of Japanese anchovy (Engraulis japonicus) in the central and southern YS move northward (Liu et al., 2020). Therefore, ENSO modifies the YS environment and impacts the ecosystems of this area.

5 Conclusions

In conclusion, we examined inter-annual changes in the BW of L. beka and A. fangsiao in the YS and used sliding window analysis and mixed-effects models to determine the impact of environmental factors on the BW of both species. We identified obvious inter-annual changes in the BW of both species. Relatively rapid environmental changes in the YS, especially SST and SOI, were reflected in the BW of both species, thereby indicating the impacts of the regional environment associated with basin-scale climate change, such as ENSO events. SST and SOI were consistently significantly correlated with BW between species, which indicates that L. beka and A. fangsiao can be appropriate indicators of environmental change in the YS. Additionally, Chl-a concentration was strongly negatively correlated with the BW of L. beka, possibly because of interspecies competition and growth rate variations. Climate change affects the catch environment, growth, and fluctuations of cephalopod populations in the YS. Knowledge of the long-term patterns of the BW of L. beka and A. fangsiao and their consistent response to environmental variables contribute to our understanding of how climate change affects cephalopod populations in the YS.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Nos. 41930534 and 41861 134037).