Interannual variability of body size and beak morphology of the squid Ommastrephes bartramii in the North Pacific Ocean in the context of climate change

Oceanic squids are typical species that sensitive to the environment change. Previous studies on oceanic squids mainly focus on its annually fluctuated abundance under the background of climate change. The changes of individual morphological sizes, including body and beak, have been observed in recent years. In this study, Ommastrephes bartramii, an important cephalopod species in North Pacific Ocean, has been selected to analyze the annual morphological changes of body and beak under different scales of climate events. Geometric morphometrics was used to analyze the shape variations of both upper beak and lower beak. The possible phenotypic plasticity of body and beak was also discussed in different aspects. Body size showed different variations among different climatic years. The size at the maturity stage in 2015 (El Niño year) was much smaller than that in other years for both sexes. The centroid size representing the beak size showed the significant difference between two types of Pacific Decadal Oscillation phases. The shapes of upper and lower beaks showed significant differences between small-scale climatic patterns in which El Niño or La Niña event occurred, displaying different growth patterns.


Introduction
Pelagic cephalopods, mostly belonging to Ommastrephidae, are a typical marine invertebrate widely distributed in the temperate and subtropical waters around the world (Chen et al., 2008;Jereb & Roper, 2010). Generally, these squids have a short lifetime less than two years (mostly one year) and are characterized by the life cycle of ''grow fast, die young'' (Pecl et al., 2004). As a result, the highly fluctuated abundance of pelagic squid has been mentioned by some researchers in recent years (Igarashi et al., 2017;Yu et al., 2019a;Alabia et al., 2020). The potential response of pelagic squid to the climate change has also been discussed (Yu et al., 2019a, b).
With the increasing number of extreme climate events, changes of feeding habitat and food composition increasingly influenced squids (Yu et al., 2019a;Portner et al., 2020). Sea surface temperature (SST) is the most frequently discussed index in previous studies (Chen, 2010;Frawley et al., 2019). The abundance of squids (always represented by catch per unit effort (CPUE) in fishery) showed significant annual fluctuations due to small-scale climate events (El Niño/La Niña) (Yu et al., 2018;Ji et al., 2020). Some large-scale climate events (e.g. Pacific Decadal Oscillation, PDO) also significantly affected marine species in some areas (Kim et al., 2018). In addition, the significant increase in SST also strongly changed the distribution pattern (Yu et al., 2019a), body size and life span of oceanic squids (Arkhipkin et al., 2015a;Takahara et al., 2017). However, the interannual variation of squid body size was rarely reported . Beak, as the main feeding organ of squid, is one of the important hard corrosion resistance structures in cephalopod (Fang et al., 2016a, b). This hard structure is not only taxon-dependent, but also stock-and sexdependent in some species (Martínez et al., 2002;Bolstad, 2006;Chen et al., 2012). The lower beak morphology or the relationship between its morphometric and mantle length is a standard criterion for the cephalopod taxonomic recognition (Jackson & McKinnon, 1996;Jackson et al., 1997;GroÈ ger et al., 2000;Lefkaditou and Bekas, 2004;Fang et al., 2015). Due to the strong relationship between the body size and beak shape, beak is a proper material for the detection of squid body change. However, the beak shape is also influenced by the changes in the ambient environment, including the habitat suitability and food availability (e.g. Dosidicus gigas (d'Orbigny, 1835), Hu et al., 2019;Ommastrephes bartramii (Lesueur, 1821), Fang et al., 2020).
Most of previous studies were focused on beak morphometrics with the traditional linear distance method (Fang et al., 2014;Liu et al., 2015). This is a simple method, although relevant measurements may increase the potential risk of bias to some degree (Bookstein, 1998;Adams et al., 2004). Moreover, the lack of geometric information also limits the application of the traditional method. Geometric morphometrics (GM) has been proposed to avoid the limitation . A large number of variables in GM can ensure that even a small shape difference (e.g. reduced morphological variability or population level) can be detected (Collyer et al., 2015). In addition, the object size is closely related to its shape (Klingenberg, 2016), and shape variation may be influenced by allometry. GM can mathematically separate shape from size and thus elucidate the level of correlation between two traits (Rohlf and Slice, 1990;Mitteroecker et al., 2013).
Ommastrephes bartramii is the most important oceanic squid widely distributed in the North Pacific Ocean (Bower and Ichii, 2005). The environmentinduced abundance fluctuation was reported (Chen, 2010;Yu et al., 2017), but the potential change in body size or beak shape is still unknown. In this study, in order to investigate the change patterns of body size and beak shape under different environment conditions, we selected squid samples from different climatic years and statistically compared the interannual body size change and beak shape difference according to GM method for the analysis of beak morphology.

Sampling and image acquisition
The samples were collected in 2009, 2010, 2015 and 2016 by the Chinese commercial jigging vessels in the main fishing ground for Ommastrephes bartramii (150°E to 160°E, 35°N to 45°N) ( Table 1). Different climate events were observed in the North Pacific Ocean in the above four years. El Niño and La Niña events were all presented by the mean values of SSTA during the sample months. The PDO index is the leading empirical orthogonal function (EOF) of monthly sea surface temperature anomalies (SSTA) over the North Pacific (poleward of 20°N) after the global average sea surface temperature has been removed (Deser et al., 2010). During a ''positive'' phase of PDO, the west Pacific becomes cooler and part of the eastern ocean warms; during a ''negative'' phase of PDO, the opposite pattern occurs (Mantua et al., 1997). In this study, large-scale climate varied from cold (negative) PDO phase in the years of 2009 to 2013 to warm (positive) PDO phase in the years of 2014 to 2016. El Niño (2009 and La Niña (2010 and 2016) also happened in the changing climatic years (Fig. 1). Therefore, the squid samples in El Niño and La Niña years from different largescaled climate events were collected and selected for morphometric analysis (Table 1).
Dorsal mantle length (ML) was measured to the nearest 1 mm after samples were captured and randomly selected on the deck. The squid samples were dissected with scissors to observe the maturity stage according to the method by Lipiński & Underhill (1995). All the upper (UB) and lower beaks (LB) were successfully dissected from the buccal mass of the samples and then stored in 75% ethyl alcohol. Since a beak was a parallel structure, we focused on one side of the beak to analyze and simulate the morphological variation. Each picture of UB and LB was photographed with Nikon D7000 camera in the same lens and a scale (ruler) was placed in each picture to enhance the morphometric standardization (Crespiabril et al., 2010;Fang et al., 2017Fang et al., , 2018. Landmarks were chosen according to the previous study, including 8 landmarks and 2 semi-landmarks for UB and 10 landmarks and 10 semi-landmarks for LB (Fang et al., 2017, Fig. 2). The combination of landmarks and semi-landmarks was used to describe the beak characteristics and better simulate the shape of beak (Neige and Dommergues, 2002;Gunz and Mitteroecker, 2013).

Geometric morphometrics
These series of landmarks were digitized by two different persons to minimize the measurement errors and the first digitized data set was chosen when the error between two measurements was less than 5% (Viscosi and Cardini, 2011). If the error was higher than 5%, a third person would be digitized again until the error between two data sets was less than 5% and then chose the first of them for the analysis (Viscosi and Cardini, 2011). A series of coordinate-based datasets were standardized in size and the landmark configuration was translated and rotated through the generalized least-squares Procrustes analysis (GPA, Rohlf and Slice, 1990). The positions of semilandmarks were optimized with Procrustes distance criterion to minimize the Procrustes Distance across samples during GPA (Bookstein, 1998). The semilandmarks slid separately, so that their positions were not influenced by the rest landmarks .  (Dryden and Mardia, 1998) and the beak shape variation associated with each principal component (PC) was graphically depicted. Multivariate patterns of beak shape variation were visualized with the scores from multivariate regressions of beak shape versus log(Centroid Size). Furthermore, in order to facilitate the climate comparison, the predicted values, which was the regression of shape on size were chosen to compare the allometric trend among factors (Adams   (Fang et al., 2017) and Nistri, 2010). Thus, the plots of the predicted values from multivariate regressions by years (in different climatic patterns) versus log(Centroid Size) were also generated (Adams and Nistri, 2010;Kelly et al., 2013). Thin-plate spline (TPS) deformation grids for the phenotypic mean shapes of O. bartramii at different scales of climatic patterns were also presented to improve the ecological interpretation of allometric changes in beak shape. All the statistical analyses were performed in R 3.5.3 (R Core Team, 2019) with the package ''geomorph'' (Adams & Otarola-Castillo, 2013).

Body size variation in different years
Mantle length showed the different changing trends among different climatic years (Fig. 3). As for females, mean ML was the largest ( The body size at the maturity stage of female squid samples in 2009 and 2010 were large than that in other two years (Fig. 3B). The individuals in 2015 were the smallest with the mean size less than 300 mm at the maturity stage. In all the male samples, the mean size was still the largest in 2009 and the smallest (292.2 mm) in 2015 and the mean sizes in 2010 and 2016 were close to 300 mm (Fig. 3B). The ML of both males and females at the maturity stage showed significant differences among different years (ANOVA, $:F = 26.91, P \ 0.0001; #: F = 23.17, P \ 0.0001).

Beak size variation in different years
The centroid size (CS) representing the beak size was calculated by GPA. The result showed that the squid beak size in cold PDO pattern (2009 and 2010) was significantly different from that of warm PDO pattern (2015 and 2016) for UB (ANOVA, F = 14.881, P = 0.0005), but the difference in LB was not significant (ANOVA, F = 0.643, P = 0.0980). The beak size in El Niño years (2009 and 2015) was in the similar small range, but more outliers (below lower quartile) were observed in warm PDO pattern (Fig. 4). The wider quartile of the beak size was depicted in cold PDO pattern with a lower median value in La Niña years (2010 and 2016) for both UB and LB (Fig. 4).

Interannual shape variation of upper beak
Sex and climatic pattern represent the covariates for the analysis of beak shape differences. As for the UB shape, sexual dimorphism was not observed based on the results of MANCOVA (Table 2). However, the UB shape difference in small climatic patterns (El Niño and La Niña events) was detected (Table 2) most probably due to the larger width of the hood observed in thin-plate spline grid configurations (Fig. S1). Allometric trends of UB growth significantly differed between small climatic patterns (Table 2). Large climatic patterns (cold PDO and warm PDO) also showed an interactive effect on small climatic pattern for UB shape ( Table 2).
The first four components (PC1, PC2, PC3 and PC4) explained 53.6% of total variation in UB shape ( Fig. 5A and B). UB shape significantly increased with the increase in body size in cold PDO pattern (Fig. 6A, left panel), but the change trend of UB shape with the increase in body size was not significant in warm PDO pattern (Fig. 6A, right panel). Predicted values showed negative trends in all the sampling years (Fig. 6B). However, the predict values of the individuals were negative/positive in cold/warm PDO pattern, respectively (Fig. 6B). The Procrustes MANOVA with sex and climatic pattern as covariates revealed that LB shape had no significant difference between sex, but the difference in LB shape between small-scale climatic pattern (El Niño and La Niña event) and large-scale climatic pattern (cold PDO and warm PDO) was significant (Table 3). The configuration result of thin-plate spline grids indicated that the lateral wall shape and width of the wing were the characteristics for distinguish the LB shapes under small-scale climatic pattern and large-scale climatic pattern, respectively (Fig. S2). Allometric growth of LB shape existed among smallscale climatic patterns as well as large-scale climatic patterns (Table 3). Large-scale climatic patterns also  had an interactive effect on small-scale climatic patterns for LB shape (Table 3).
The first four components (PC1, PC2, PC3 and PC4) explained 55.3% of total variation in LB shape ( Fig. 7A and B). The LB shape varied with the  (Fig. 8A). Predicted values showed positive trends in all the sampling years (Fig. 8B). Most of the individuals' predicted values in cold PDO pattern were negative (Fig. 8B, left panel) and some of the predicted values in warm PDO pattern (most individuals in 2016 and half of individuals in 2015) were positive (Fig. 8B, right panel).

Body size differences
As a typical phenotypic plasticity species, squids always tend to change their morphology to adapt the changing environment (Arkhipkin et al., 2015a;Jones et al., 2019). In the selected years of this study, the high proportion of small ML groups (\ 240 mm) of O. bartramii occurred in 2010 (more than 50% of small groups, Fig. 3) and 2015 (nearly 40% of small groups) and no male sample [ 360 mm (Fig. 3). In El-Niño/ La Niña years, the small body size could be a consequence of extreme high/low environmental temperatures for other oceanic squids (D. gigas, Keyl et al., 2011), as observed in this study. A higher ambient temperature led to the shrinkage of squid body size, as reported in other species (D. gigas: Arkhipkin et al., 2015a; Todarodes pacificus (Steenstrup, 1880): Takahara et al., 2017). The annual body size variation of D. gigas, which is a sibling species with O. bartramii and has similar biological features (Arkhipkin et al., 2015b), had been discussed in previous studies Keyl et al., 2008Keyl et al., , 2011Arkhipkin et al., 2015a). Environmental conditions and food availability are the two main causes for this interannual change (Keyl et al., , 2011. Unfavorable temperature might restrain the metabolism of squids, so this phenomenon might also promote body size miniaturization, just like the size at the maturity stage in 2015 in this study (Fig. 3B).
The annual difference in body size of female samples was much bigger than that of male samples (Fig. 3). Body size of females tended to be more plastic than that of males (Fig. 3). Most of previous studies also discussed the variation of female body size in other squid species (e.g. D. gigas, Argüelles et al., 2008;Keyl et al., 2008Keyl et al., , 2011Arkhipkin et al., 2015a). In general, female squids grow faster than males due to their high-rate feeding, thus leading to the faster metabolism and gonadal growth (Markaida et al., 2004;Fang et al., 2016b). This was shown for O. bartramii (Fang et al., 2016b) and D. gigas (Markaida et al., 2004), but also in other oegopsid species, such as Illex illecebrosus (Lesueur, 1821) and Illex argentinus (Castellanos, 1960) (O'Dor 1998. This is a kind of survival strategy to maintain the population in a certain scale (O'Dor 1998).

GM-based beak shape variation among different years
Geometric morphometrics (GM) is a useful method to examine the squid beak shape variation at the stock level (Fang et al., 2017) or in different years (this study), but such variations are too small to be detected. The unusual small centroid size outliers in 2015 ns Not significant *P significant at a = 0.05 (Fig. 4), which were coincident with the ML distribution (Fig. 4), indicated that these small beak shapes were closely related to the slower growth in the beak.
Different scales of climatological and environmental variability, including ENSO events (Chen et al., 2008) and different PDO phases (Yu et al., 2017), may  (Villanueva et al, 2017). This variation of feeding habits, including prey item, may also act as one of the main factors influencing the squid growth, as well as the beak growth. The results of statistical analyses in this study suggested that the beak shape of winter-  (Tables 2 and 3). The change of feeding habits (i.e. prey type) for squid might be one of the reasons for the variation of beak shape, as reported in other cephalopod species (Matias et al., 2019). The annual beak shape variation was difficult to be identified only from the reconstructed TPS grids (Figs. S1 and S2) because this annual environmental change might not substantially affect the beak shape. However, GM-based statistical analysis results indicated that the annual difference was significant.
Relationship between oceanic squid morphological variation and climate change In the last several decades, climate change events frequently occurred and induced the landing variation of squids (Igarashi et al., 2017). Actually, the landing fluctuation and the morphological change had been discussed for other squid species in previous studies Keyl et al., 2008Keyl et al., , 2011. The body size change of O. bartramii in El Niño year was similar to that of other Ommastrephidae species, thus it is important to understand the beak morphological variation, which is closely related to the feeding behavior (Hu et al., 2019). Beak shape variation was consistent with the mantle length change in this study, so it is reasonable to emphasize the impact of feeding behavior on the squid body size. Indeed, the morphological size would change under the background of climate change. However, the growth pattern change was seldom reported. Here, in this study, we found that the growth pattern of beak shape was changed in different climate events, regardless of UB and LB (Tables 2 and 3). Oceanic squid is a typical opportunist species and has a large feeding spectrum. A squid can utilize surrounding foods as possible to keep the fast growth rate, even in the circumstance of food deficiency (Watanabe et al., 2004(Watanabe et al., , 2008. In this study, O. bartramii both changed its body size and beak shape. In addition, the body size had a strong resilience and could be largely changed within a short time. This phenomenon was also observed in the beak morphology (Figs. 4 and 5). Thus, the phenotypic plasticity of O. bartramii is a kind of response to the environmental change and can be regarded as a survival strategy.