Trophic ecology of Humboldt squid, Dosidicus gigas, in conjunction with body size and climatic variability in the Gulf of California, Mexico

Over the past two decades, the Gulf of California (GOC) has experienced three strong El Niño events (1997–1998, 2009–2010, and 2015–2016), each of which was followed by a drastic reduction in mantle length of mature Humboldt squid, Dosidicus gigas (from >60 cm to <20 cm). However, it is unclear how the oceanographic changes associated with strong El Niño events affected the midwater organisms on which D. gigas feed, limiting our ability to assess the relative importance of temperature and food availability in the phenotypic response of D. gigas to environmental variability. We quantified the diet of D. gigas in the GOC before, during, and following the past three El Niño events and found that although its diet varied little across a large range of body sizes (8–85 cm), significant and predictable diet variability was observed with respect to sea surface temperature and chlorophyll‐a concentration. Consumption of large numbers of relatively small, high calorie prey in both relatively cool (anchovies) and relatively warm, productive conditions (myctophids) is likely necessary to support growth to large body sizes before maturation. When warm, unproductive conditions prevailed in the GOC, only small squid were present and had diets dominated by euphausiids and pteropods, prey with relatively low caloric value. Using a time series of diet data, this work provides unique insights into the response of a midwater forage community to oceanographic variability and the effects of environmental variability on the trophic ecology of an oceanic predator.


Abstract
Over the past two decades, the Gulf of California (GOC) has experienced three strong El Niño events (1997-1998, 2009-2010, and 2015-2016), each of which was followed by a drastic reduction in mantle length of mature Humboldt squid, Dosidicus gigas (from >60 cm to <20 cm). However, it is unclear how the oceanographic changes associated with strong El Niño events affected the midwater organisms on which D. gigas feed, limiting our ability to assess the relative importance of temperature and food availability in the phenotypic response of D. gigas to environmental variability. We quantified the diet of D. gigas in the GOC before, during, and following the past three El Niño events and found that although its diet varied little across a large range of body sizes (8-85 cm), significant and predictable diet variability was observed with respect to sea surface temperature and chlorophyll-a concentration. Consumption of large numbers of relatively small, high calorie prey in both relatively cool (anchovies) and relatively warm, productive conditions (myctophids) is likely necessary to support growth to large body sizes before maturation. When warm, unproductive conditions prevailed in the GOC, only small squid were present and had diets dominated by euphausiids and pteropods, prey with relatively low caloric value. Using a time series of diet data, this work provides unique insights into the response of a midwater forage community to oceanographic variability and the effects of environmental variability on the trophic ecology of an oceanic predator.
Warming of the global ocean is affecting oceanic ecosystem dynamics by altering the distribution of marine habitats (Doney et al. 2011;Hazen et al. 2012), shifting the sizestructure of marine communities (Polovina and Woodworth-Jefcoats 2013;Reuman et al. 2014), and altering predator-prey dynamics (Green and Côté 2014). However, these ecosystemwide effects are rooted at the organismal level. Warming increases metabolic demand, variability in growth rate, timing of maturation, and fecundity (Houde 1989;Forsythe et al. 2001;Gillooly et al. 2001;Pörtner et al. 2001;Pörtner and Peck 2010). Changes in size of both predators and prey may alter the composition and size structure of accessible prey communities through size-based constraints on predation (e.g., gape limitation) and behavior (e.g., vertical distribution) (Scharf et al. 2000;Heupel et al. 2014). Climate-driven changes in prey selectivity as well as prey availability thus both impact oceanic ecosystems. Variability in prey selectivity during ontogeny and growth, as well as under different environmental conditions, is often assessed through diet analysis. Understanding the dynamic relationships between oceanic predators and their prey is particularly important for the sustainable use of marine resources in a changing climate (Pikitch et al. 2004).
Dosidicus gigas, an ommastrephid squid commonly known as the Humboldt or jumbo flying squid, is an oceanic predator that supports the largest single-species invertebrate fishery in the world (FAO 2018). D. gigas is also an important food source for many economically valuable fishes such as sharks and billfishes (Preti et al. 2012;Rosas-Luis et al. 2016), and for many marine mammals, including sperm whale and porpoises (Perrin et al. 1973;Clarke and Paliza 2001). Throughout its range in the eastern Pacific Ocean, D. gigas consumes vertically migrating, midtrophic organisms that occupy mesopelagic depths (200-1000 m) during the day and migrate to near-surface waters at night, particularly myctophid and phosichthyid fishes (Markaida and Sosa-Nishizaki 2003;Alegre et al. 2014).
Across its range, D. gigas has relatively low genetic diversity (Sandoval-Castellanos and Uribe-Alcocer 2007;Sandoval-Castellanos et al, 2010;Staaf et al. 2010) but exhibits large phenotypic variation in mature size (13.5-120 cm mantle length, <0.1 to >60 kg), most obviously with respect to latitude (Nigmatullin et al. 2001). Individuals living in subtropical or temperate, productive habitats like the California Current System and Humboldt Current System grow for up to 2 years to large sizes (>60 cm mantle length) before maturing and entering their terminal spawning period (Arkhipkin et al. 2015). Conversely, D. gigas in tropical, less productive regions, such as the Costa Rica Dome, mature and reproduce in less than 1 yr at much smaller sizes (<20-40 cm mantle length) (Chen et al. 2014;Gong et al. 2018).
Phenotypic plasticity in this species also occurs between cohorts within geographically defined populations, particularly in the Humboldt Current System, where oceanographic circulation and productivity patterns are complex (Argüelles et al. 2001;Chen et al. 2011). In the Guaymas Basin of the Gulf of California (GOC), a comparatively small body of water, uniformly large squid historically supported a substantial commercial fishery starting in the mid-1990s. Strong El Niño events in 1997-1998, 2009-2010 were all followed by a reduction in the mature size of D. gigas in the next generation from~60 cm to~20 cm mantle length (Hoving et al. 2013;Robinson et al. 2016;Gilly, unpubl.). In each case, the absence of large squid led to fishery collapse. Although both mature size and commercial landings recovered after El Niño 1997-1998, there has been essentially no recovery of either since El Niño 2009-2010 (Robinson et al. 2016;Gilly unpubl.), suggesting that the intervening conditions between El Niño events are critical for the recovery of squid size in the GOC. Ocean warming and low primary productivity, signatures of a strong El Niño, are clearly associated with the change in size of D. gigas in the GOC to the small, equatorial phenotype (Hoving et al. 2013). Persistence of these environmental conditions since El Niño 2009-2010 is thus also associated with the lack of recovery in size at maturity during this period (Robinson et al. 2016;Frawley et al. 2019;Gilly et al., unpubl.). This history presents an opportunity to examine how climatic variability can influence the diet and inferred prey field of D. gigas, a factor that is undoubtedly relevant to the persistence of small mature D. gigas in this area.
Although there have been several dietary studies of D. gigas in the GOC over the past 20 years (Markaida and Sosa-Nishizaki 2003;Markaida 2006;Markaida et al. 2008;Hoving et al. 2013), none has explicitly assessed the relationship between environmental conditions and diet variability nor any patterns of diet variability across the large size range of adult squid observed. The goals of this study were (1) to quantify the diet of D. gigas in the GOC during the 2015-2016 El Niño, (2) to assess whether its trophic ecology varies with respect to body size and environmental variability by comparing our data to historic diet studies from the same area, and (3) to explore whether D. gigas diet composition can be used to infer prey availability in the GOC.

Data and methods
Squid sampling and diet analysis D. gigas was captured with jigs (4-15 cm long) near the surface using rod-and-reel and handlines from small fishing vessels near Santa Rosalía, BCS, Mexico, in July 2015 (n = 96, from 5 sampling events) and throughout the central GOC from the R/V El Puma during oceanographic cruises in June of 2016 (n = 196, from 14 sampling events) and 2017 (n = 100, from 15 sampling events) (Fig. 1A). Dorsal mantle length, sex, and maturity state (Lipi nski and Underhill 1995) were recorded for all individuals (summarized in Fig. S1). Whole stomachs were excised and frozen at −20 C. Stomachs were later thawed and weighed to the nearest 0.01 g. Contents were removed from the stomach and rinsed over a 250 μm Nitex mesh filter. Empty, rinsed stomachs were reweighed to determine the wet weight of the contents. Stomach contents were transferred to an 8.5 cm petri dish and diluted in tap water (to a final volume of 20-30 mL) for sorting and quantification. If the total content weighed more than 1 g, it was subdivided into multiple petri dishes (with <1 g of contents per dish) to facilitate thorough sorting.
Each sample was examined with a binocular dissecting microscope (8×-40×) under both dark-field and bright-field illumination. Prey items much smaller than the expected prey size for squid of the sizes sampled (<0.01 g, Fernandez Araoz 1991; Barnes et al. 2010;Romero-Romero et al. 2016) and commonly consumed by myctophids, including ostracods, copepods, and amphipods (Dalpadado and Gjøsaeter 1988;Van Noord et al. 2016), were considered to have been prey of consumed fishes and excluded from the analyses.
Fish otoliths, cephalopod beaks, pteropod shells, and crustacean parts were identified to the highest possible taxonomic resolution following Markaida and Sosa-Nishizaki (2003). The frequency of occurrence for each prey item was calculated as the percentage of squid stomachs containing that prey item (Markaida and Sosa-Nishizaki 2003). Fish and cephalopod prey were enumerated as the maximum number of right or left fish otoliths, and of upper or lower cephalopod beaks, respectively. Pteropods were enumerated as the number of shell apices, but we were not always able to accurately count crustaceans. For the purpose of estimating their numerical abundance in the diet, crustaceans were recorded as present (n = 1), or absent (n = 0). Although this protocol likely underestimates the numerical importance of crustaceans, it provides a conservative estimate of crustacean consumption. Squid heavily masticate their food, making it difficult to accurately separate and weigh all tissues by prey type, eliminating the possibility of directly measuring the gravimetric contribution of prey types to their diet. Numerical composition is a more conservative estimation of diet composition for squid and is a more useful metric for examining prey encounter rates as a proxy for prey availability.
To address the possibility that changes in trophic ecology may be reflected as changes in prey size, body length was estimated for the two most common fish and most common cephalopod prey species. Fish sagittal otoliths and rostra of cephalopod beaks were measured to the nearest 0.01 mm using an ocular micrometer. Standard lengths for the prey fishes Benthosema panamense and Triphoturus mexicanus were calculated based on published regressions of otolith lengths (Markaida and Sosa-Nishizaki 2003;Markaida et al. 2008). Mantle lengths of Pterygioteuthis hoylei were calculated based on a published regression of beak rostral lengths (Wolff 1984, reported as Pterygioteuthis giardi). Only minimally digested sagittal otoliths and cephalopod beaks with intact rostra were used in our analyses to reduce the chances of underestimating prey size.

Combining new diet data with previous studies
Data collected during 2015-2017 were combined with four, published D. gigas diet datasets collected from our study area during 1996-2011 (Markaida and Sosa-Nishizaki 2003, Markaida 2006, Markaida et al. 2008, Hoving et al. 2013 ,  Table S1). Squid collection and diet analysis protocols from these studies were similar to our own, but their taxonomic resolution and prey treatment were adjusted to match the 2015-2017 dataset to enable diet comparison across the entire time series. All pteropods and crustaceans identified to a higher taxonomic resolution in prior datasets were aggregated into "pteropod" and "crustacean" groups to account for lack of further identification of these groups during 2015-2017. For stomachs containing more than one crustacean prey (n = 72, 8.8% of stomachs containing prey from previous datasets), numerical abundance was converted to the presence/absence as described above. Items considered to be prey of consumed myctophids, rather than of D. gigas, were excluded from previous datasets prior to analysis. Pieces of flesh identified as recently ingested D. gigas based on tissue condition and thickness (~1% of total numerical composition) were considered fishing-associated cannibalism and excluded from our analyses.
Standard lengths for B. panamense and T. mexicanus, and mantle lengths for P. hoylei were estimated using the protocols described above and reported in Markaida and Sosa-Nishizaki (2003) and Markaida (2006), while those from Hoving et al. (2013) are estimated and reported here for the first time (summary of prey length data and their origins are presented in Table S2).

Environmental conditions in the GOC
Sea surface temperature (SST), surface chlorophyll a concentration (CHL, as a proxy for phytoplankton biomass), and NOAA's Multivariate ENSO Index (MEI, Wolter and Timlin 1993) were used to quantify environmental conditions within our study area in the GOC. SST is considered here as a proxy for conditions experienced by D. gigas and its prey, many of which migrate into near-surface waters at night. Satellite observation-based SST and CHL products were accessed from the E.U. Copernicus Marine Service Information (http://marine. copernicus.eu). Monthly averages of SST at 0.25 × 0.25 resolution covering the entire diet time series were obtained from the Global ARMOR3D L4 dataset through the MULTIOBS_GLO_ PHY_REP_015_002 product. Monthly averages of CHL at 4 km × 4 km resolution were obtained from the OCEANCOLOUR_GLO_ CHL_L4_REP_OBSERVATIONS_009_082 product for September 1997-December 2016 and from the OCEANCOLOUR_GLO_ CHL_L4_NRT_OBSERVATIONS_009_083 product for January-August 2017. These data products were selected to maximize the temporal coverage of environmental data using the minimum number of products. Bimonthly MEI values were obtained from the NOAA Earth Systems Research Laboratory (https://www.esrl. noaa.gov/psd/enso/mei/table.html) and were designated as the MEI value for second month of the index (e.g., the "DECJAN" index was treated as January, etc.).
To minimize the effects of variable spatial extent of squid sampling between years on diet similarity, we split our study area into five regions with respect to oceanography and biogeography (northern, Baja shelf, central pelagic, Sonora shelf, and southern; Fig. 1B). The northern region, separated from the rest of our study area by a relatively shallow sill of several hundred meters depth (Alvarez-Borrego and Lara-Lara 1991; Gilly et al. 2012), experiences tidally driven upwelling (Roden 1964) and supports a distinct community of micronektonic fishes (Robison 1972). The southern region includes the Carmen Basin and approximates biogeographic boundaries proposed by Gilbert and Allen (1943) and Round (1967). The central gulf, mostly encompassing the Guaymas Basin was split into three regions using the 1000 m isobath to reflect seasonal variability in upwelling (off the mainland coast in winter and off the coast of the Baja California peninsula in summer, Roden 1964), distinct micronektonic fish and euphausiid communities (Robison 1972;Lavaniegos et al. 1989), and potential variability in prey assemblages between shelf and pelagic habitats as seen in other ecosystems (Reid et al. 1991;Alegre et al. 2014).
Monthly environmental data were grouped to obtain regional averages of SST and CHL. Regional averages at monthly lags from 0 to 18 months were also included in the analyses to allow for the possibility of a delay in response time of prey communities to variability in surface habitat temperature and phytoplankton biomass. The average maximum life span of common prey of D. gigas was estimated to be approximately 18 months based on studies of myctophids and euphausiids from other subtropical regions (Gartner 1991;Harvey et al. 2009) and was used to define the maximum lag included in our analyses.

Classification and regression tree (CART) analysis and diet similarity
Diet composition was quantified as the numerical proportion of prey items within an individual stomach. Prey types that accounted for more than 1% of the total numerical proportional abundance were included in the analyses (after Olson et al. 2014;Portner et al. 2017), and prey types that contributed less than 1% to the total proportional abundance were grouped together by a broad taxonomic identifier ("fish," "cephalopod") and treated together as an "Other" prey group. Because we were interested in testing the effects of variability of specific habitat characteristics on squid diet, we excluded capture year as an explanatory variable from our analyses, which aggregates habitat variability into a single, confounded variable. Capture months were grouped by "season" (warm = June-October, cool = November-May) based on GOC-wide trends in SST and historic squid migration patterns in the Guaymas Basin (Markaida et al. 2005;Boyer et al. 2013).
To examine the effects of capture region, season, SST, CHL, MEI, mantle length, sex, and maturity state on diet composition of D. gigas, we used classification and regression tree (CART) analysis in the diet package (Kuhnert and Duffy 2013) in R (R Core Team 2015). For a more thorough discussion of this method and its application to diet studies, see Breiman et al. (1984) and Kuhnert et al. (2012). 10-fold cross validation was performed on the full tree, which was then pruned to a size that minimized the cross-validated prediction error rate. We selected the least-complex, pruned tree within one standard error (SE) of the tree with the minimum cross-validated error rate following the "1-SE" rule (Breiman et al. 1984). The importance of each variable for explaining diet variability was determined as described in Kuhnert et al. (2012) and are reported relative to the explanatory variable that most accurately predicted the splits at all nodes throughout the tree. Diet diversity is reported as the Gini Index of Diversity (D) for each terminal node of our selected, pruned regression tree (Breiman et al. 1984).
Variables describing both past and present environmental conditions were selected to account for differences between short-and long-term effects of habitat conditions on prey availability (e.g., recruitment vs. habitat selection). Lags of environmental variables used to describe past conditions were selected for the final CART analysis by building separate CARTs using all variables at a single lag and comparing classification performance (quantified as classification error rate) between the trees (Table S3). For each environmental variable, the lags producing the two best-performing trees where that variable had the highest relative importance were considered for inclusion in the final CART. For SST and MEI, the best-performing trees were clustered within three consecutive months of consistently high-performing trees with similar variable importance (Table S3). In these cases, all three lags were considered for inclusion in the final CART regardless of performance-rank to capture the consistent explanatory power of SST and MEI during these periods. SST and CHL at 0-and 1-month lags were considered as present conditions. This method resulted in two possible pairs of variables defining present environmental conditions and 24 possible combinations defining past environmental conditions. CARTs using all possible combinations of single-lag present conditions and mixed-lag past conditions were built and compared to determine the combination of variables that explained the most variance within our diet dataset. Lags are reported here as the two-digit number of months prior to the month of capture following the environmental variable (e.g., SST.01 is the mean SST from the region of capture 1 month prior to the month of capture).
Diet compositions summarized and discussed are the mean cross-validated proportions at each terminal node. Traditional bootstrap aggregation (n = 500) was performed on our final tree in diet. Though each squid was captured individually, squid from the same sampling event (defined as a single, stationary fishing session) may not be completely independent, as multiple individuals were presumably captured from the same group during single fishing events. We examined whether utilizing multiple individuals from the same sampling event introduced pseudo-replication into our analyses by constructing an experimental variogram of spatial bootstrap prediction residuals (Hunsicker et al. 2012;Kuhnert et al. 2012).
Diet similarity between samples within nodes identified by the CART analysis was quantified using the Morisita-Horn Dissimilarity Index (Ĉ MH , Horn 1966) in Primer v7 (Clarke and Gorley 2015).Ĉ MH is an appropriate metric of compositional similarity between groups with different sample sizes for data reported as proportional abundance (Krebs 1999). We performed permutational multivariate ANOVA (PERMANOVA; Anderson et al. 2008) on the diet dissimilarity matrix in PRIMER v7 (Clarke and Gorley 2015) to examine whether specimens in different nodes had significantly different diet compositions. To account for differences in sample size between nodes, we used Type III sum of squares (SS) for the PERMANOVA (Anderson et al. 2008).

Prey length
The relationship between prey length, D. gigas mantle length, capture region, and environmental conditions was examined using multiple linear regression models in R (R Core Team 2015). To minimize the number of interaction terms and prevent overfitting, we used a single variable to represent environmental conditions for each model. Models with different environmental variables and lags were compared using Akaike's Information Criterion, and the most parsimonious model was selected for each prey species. Only prey from D. gigas collected in the warm season (June-October) were included in our analyses of prey lengths to minimize potential effects of season on population size-structure (Smoker and Pearcy 1970;Gjøsaeter 1987). Analyses were also restricted to individuals from the central GOC (Regions B, C, D, Fig. 1B), due to low sample size of measured prey in the southern and northern regions. For all models, ANCOVA was performed with Type III SS using the car package (Fox and Weisberg 2011), and significant interactions between explanatory variables of the linear models were visualized with partial residual plots using the visreg package (Breheny and Burchett 2017).

Diet description for squid collected 2015-2017
Of the 392 stomachs from D. gigas captured during 2015-2017, 9.95% were empty, 14.8% contained D. gigas soft tissue, and 14.8% contained inorganic items (mostly small plastic fibers of various colors), although an unknown fraction of these items may have been initially ingested by prey. For the 353 stomachs containing material considered to be prey, diet composition was dominated by the myctophid B. panamense, crustaceans, and pteropods. We observed 1007 fishes, 149 cephalopods, 246 pteropods, and a minimum of 227 crustaceans. Of the 1007 fish prey, 211 were treated as "unidentified," because they were represented by skin, non-otolith bones, unidentifiable fragments of otoliths, or relatively small otoliths that could not be identified beyond "non-myctophid." Two prey fishes, referred to as "Species A" and "Species B," were associated with distinct otoliths and observed repeatedly but were not successfully identified (photographs given in Fig. S2). A summary of the abundance and frequency of occurrence for all identified prey items is given in Table 1.

Environmental variables are important determinants of diet variability: CART analysis
CART analysis was performed on diet data from 1387 D. gigas specimens collected during 1996-2017 and spanning 8-85 cm mantle length (Fig. 2). We selected the most parsimonious, pruned tree within one SE of the tree with the lowest classification error rate. This tree, based on the 1169 stomachs that contained prey, had seven splits and a crossvalidated error rate of 0.267 (SE = 0.019, R 2 = 0.146, Fig. 3A). Ten-fold cross validated predictions of mean proportional diet composition for each node are given in Fig. 4. Spatially bootstrapped predictions (n = 500) for the proportional abundance of each prey group within each terminal node are given in Fig. S3. Environmental variables that best explained diet variability in our dataset were CHL.01 (Rank = 0.42) and SST.01 (Rank = 0.79), describing conditions at the time of squid capture, and CHL.06 (Rank = 0.41), SST.11 (Rank = 0.87), and MEI.13 (Rank = 1.00) describing past conditions (Fig. 3B). When included in the CART analysis, season, sex, and maturity state explained no variability in our diet dataset (Ranks = 0) and were therefore dropped from our analyses. An experimental variogram of the spatial bootstrap prediction residuals exhibited a very small nugget effect and demonstrated a minimal effect of distance on variance (Fig. S4), suggesting that spatial autocorrelation among individuals from the same sampling event is negligible and that traditional bootstrap methods sufficiently estimate error of diet predictions (Dale and Fortin 2014). The number of stomachs and mean diet diversity for each node are reported in Table 2.
The deepest split in our tree separated D. gigas specimens based on ENSO conditions from the previous year (Fig. 3A). When the MEI 1 year prior to capture of the squid was negative (MEI.13 < −0.09), conditions were relatively cool and productive with all sampled squid being larger than 30 cm mantle length and having a diet dominated by B. panamense (Node 2, MEI.13 < −0.09, n = 250, D = 0.36) (Fig. 4). If MEI 1 year prior to capture was positive (MEI.13 ≥ −0.09), a state generally associated with relatively warm and unproductive conditions, local environmental conditions at the time of capture (SST.01, Rank = 0.79), mantle length (Rank = 0.62), and capture region (Rank = 0.22) were important determinants of diet composition (Fig. 3A). Thus, if MEI.13 was positive, squid of different size classes had significantly different diets under both cool and warm environmental conditions (SST.01). The size at which the dietary distinction occurred was similar (~25 cm mantle length) for both cool (27.3 cm; SST.01 < 21.04 C, Node 13, n = 142) and warm conditions (21.55 cm; SST.01 ≥ 21.04 C; Nodes 28 and 29, n = 394).
1.4 mg/m 3 , n = 331, D = 0.81), whereas the diet of smaller squid was dominated by crustaceans (Node 61, CHL.01 < 0.94 mg/m 3 , n = 24, D = 0.57). Small squid collected in southern, northern, and central pelagic regions (i.e., not Baja or Sonora shelf, regions B and D in Fig. 1B, respectively) had relatively diverse diets composed of fishes, cephalopods, and pteropods irrespective of local productivity (Node 31, n = 165, CHL.01 = 0.33 < CHL.01 < 4.04 mg/m −3 . Diet similarity D. gigas in different nodes of our regression tree had significantly different diets (PERMANOVA, p = 0.001-0.021, pseudot = 1.747-15.005; Table 2), but squid collected under similar environmental conditions had more similar diets than those collected under different environmental conditions, regardless of squid size. In particular, diets of squid captured when MEI.13 was negative (Node 2), characterized by relatively cool, productive conditions, were comparable to those of squid captured when MEI.13 was positive, but present conditions were relatively productive (Nodes 28 and 60). These groups were the most similar across our dataset (Ĉ MH Node 28: Node 2 = 0.54, C MH Node 60: Node 2 = 0.64,Ĉ MH Node 60: Node 28 = 0.68), yet they encompassed the largest range of squid mantle length (~22 -85 cm). When MEI.13 was positive, and present conditions were relatively warm and less productive, smaller squid had relatively similar diets in all capture areas (Ĉ MH Node 61: Node 31 = 0.65). Larger squid collected when relatively cool conditions followed MEI.13 positive conditions had a diet dominated by anchovies  (Engraulis mordax), a situation that was most distinct from all other nodes (Node 13: D = 0.47,Ĉ MH = 0.86-0.92). Larger squid collected during relatively warm and less productive conditions also had a diet composition that was quite distinct from all other nodes, but this difference was due to the mixed consumption of many common prey species, rather than to an unique prey composition (Node 29: D = 0.81,Ĉ MH = 0.81-0.92). Proportional numerical abundance and frequency of occurrence for all prey groups are reported for all terminal nodes in Table S4.

Prey size
We measured a total of 957 otoliths from B. panamense, 356 from T. mexicanus, and 517 beaks from P. hoylei collected from 1996 to 2017 between June and October. Over 95% of all measured prey were collected from the Baja, Sonora, or pelagic regions of central GOC and were included in our analyses. For all prey species measured, selected multiple linear models and ANCOVAs are given in Table 3. The number of prey measured per year is given in Table S2. A summary of all multiple linear models tested is given in Table S5. Mean prey size generally increased with increasing D. gigas mantle length, though there were specific prey-type exceptions to this trend with respect to environmental conditions and capture region. The range of prey sizes consumed for each prey species measured remained fairly constant, regardless of squid size or environmental conditions. A positive relationship between B. panamense standard length and D. gigas mantle length was found for positive MEI.01 conditions, with the slope increasing with larger values of MEI.01 (Fig. 5A). Conversely, the size of B. panamense from squid captured during ENSO-negative conditions had a negative correlation with D. gigas mantle length (Fig. 5A). However, data for MEI.01 = 1.5 were limited to small D. gigas, whereas data for MEI.01 = −05 and + 0.5 span a much larger range of mantle lengths. A positive slope for the relationship between B. panamense and D. gigas mantle length was observed for squid from Baja and Sonora shelf regions, but this relationship had a negative slope for squid from the centralpelagic region of the Guaymas Basin (Fig. S5A).
We found a weakly negative but significant relationship between standard length of T. mexicanus and ENSO conditions with a 14-month lag (MEI.14; Fig. 5B) and a positive relationship between these parameters for squid collected on the Baja shelf and central-pelagic region (Fig. S5B). This relationship was not observed for the Sonora shelf (Fig. S5B). The interaction effect of ENSO conditions and D. gigas mantle length on T. mexicanus standard length was significant but very weak and was not considered further (Table 3).
P. hoylei mantle length showed a complex relationship with D. gigas mantle length that depended on both capture region and ENSO conditions 4 months prior to capture (MEI.04; Table 3). In the central-pelagic region, a positive relationship between P. hoylei and D. gigas size was found for all ENSO condition, with the slope increasing with increasing values of  Table 2. Average pairwise Morisita-Horn dissimilarity (0 = completely similar, 1 = completely dissimilar) for individuals between nodes below diagonal, and pairwise pseudo t-statistics between nodes from the PERMANOVA above diagonal (all p ≤ 0.001 unless indicated with a "*," for which 0.001 < p < 0.05). The number of specimens (n) and the Gini Index of Diversity (D) for each node are also given. MEI.04 (Fig. 5C). For both shelf regions of Baja and Sonora, the corresponding slope changed sign between negative and positive values of MEI.04, then increased with increasing MEI-but the slopes of the two regions were always opposite in sign (Fig. 5C). Again, the size range of D. gigas sampled for MEI.04 = 1.5 was very limited.

Discussion
D. gigas exhibits significant variability in diet with respect to capture location, body size, and environmental conditions in the GOC (both at and prior to the time of capture). Differences in diet composition were found between shelf and pelagic habitats in the central GOC, as well between these areas and biogeographically distinct regions north and south of the Guaymas Basin (Robison 1972;Brinton et al. 1986;Avendaño-Ibarra et al. 2013). We sampled a large range of body lengths (<10 to 85 cm mantle length) across more than two orders of magnitude in body weight (<0.1 to 10 kg) and consistently observed a single dietary shift with respect to body size at~25 cm mantle length for squid captured when MEI was positive 1 year earlier (Fig. 3A). If MEI was negative 1 year before capture, the diet for all squid was dominated by B. panamense (Node 2, Fig. 4), though no small squid were sampled under these conditions that seem to favor larger individuals. In general, the majority of diet variability in our dataset could be explained by environmental conditions approximately a year prior to capture and at the time of capture. These findings suggest that prey availability depends on both past and present environmental conditions and is a strong determinant of diet variability throughout the life of an individual squid. Maturity state and sex of individual squid explained no variability in our diet dataset and were excluded by our analyses, suggesting that the changes in prey selectivity are driven wholly by changes in body size, independent of sex or maturity state for the squid sampled.

Size-structured shift in prey selectivity
Despite the wide range of environmental conditions experienced by squid in our dataset, there was consistent diet variability between "large" and "small" (>~25 cm and <~25 cm mantle length, respectively) D. gigas, regardless of environmental conditions. Splits in our tree that partition squid by size have individuals from the same years on both side of the split (for most years), though we sampled few squid <25 cm mantle length during years when the size of mature squid was relatively large. Squid in a given area at a given time have access to the same prey; thus, diet variability between size classes is best explained by differences in prey selectivity with respect to body size. This shift in selectivity likely reflects a combination of variability in metabolic demand, swimming ability, and size or morphology of the feeding apparatus (arms, tentacles, and mouth) with increasing body size (Rodhouse and Nigmatullin 1996;Phillips et al. 2003).
Mass-specific metabolic rate of D. gigas decreases with increasing body size, perhaps allowing for an increased allocation of energy to reproduction ). Variability Table 3. Summary of selected multiple linear regression models for each prey species. The total number of prey individuals measured for each species (n) is given as well as the number from this study (2015)(2016)(2017) and from Hoving et al. (2013Hoving et al. ( ) (2010Hoving et al. ( -2011 that were previously unreported. The formula, F-statistics, adjusted R 2 , and p value are given for each model. SS, F-statistics, and p values from ANCOVA are also given for each model. Terms in bold face for each model are visualized in Fig. 5. Benthosema panamense n (2010)(2011)(2015)(2016)(2017) 917 ( in nutritional requirements and consumption rates with respect to body size in D. gigas are not well studied. Ehrhardt (1991) estimated consumption rate of D. gigas to be 13.4% body weight per day, but the size of the squid studied was not indicated nor were the caloric values of different prey items considered. We observed that large squid ate more fishes and cephalopods, and fewer crustaceans than do small squid, representing a shift to more calorie-dense prey with increasing body size. Euphausiids typically have the lowest energy density per unit weight, cephalopods have an intermediate energy density, and fishes (myctophids and anchovies) have the highest energy density (Benoit-Bird 2004;Glaser 2010;Sinclair et al. 2015). Thus, even if consumption rate remains constant across a large range of body sizes, increased consumption of fish with increasing size would increase the total calories ingested.
There are few relevant data concerning swimming and feeding mechanisms for individual D. gigas, but both likely affect prey selectivity. Acoustic observations of free-swimming D. gigas in the GOC recorded an approximate doubling of maximum velocity between squid~25 cm and 80 cm mantle length (8-20 m/s, Benoit-Bird and Gilly 2012). These velocities greatly exceed the maximum sustained and burst speeds of common crustacean (0.06-1 m/s; O'Brien 1987, Benoit-Bird 2004), cephalopod (0.4-1.4 m/s; Young et al. 1982, O'Dor andWebber 1991), and fish prey (0.2-3 m/s; Sambilay 1990, Dypvik et al. 2012. Although swimming velocities directly recorded from large, tagged individuals of D. gigas are lower (2.5-10 m/s; Gilly et al. 2012 and Gilly et al., unpubl.), it is likely that maximum swimming speed does not constrain feeding on most common prey items observed for large or small size classes of D. gigas.
Feeding behaviors exhibited by large D. gigas observed from remotely operate vehicles appear to minimize energy expenditure for prey capture and include tentacle striking and "cage feeding," where all eight arms are spread to engulf smaller prey (Schlining and Stout 2006;Young and Vecchione 2013;Trueblood and Seibel 2014). Presumably, these behaviors also occur in the absence of artificial lighting from underwater vehicles and are similar in small squid. For any "ambush"style feeding interactions, the longer arms and tentacles of larger squid may increase the distance between predator and prey at which a successful feeding event could occur, minimizing the chances of eliciting an escape response by the prey. In both cases, squid approach the prey slowly and fast swimming is not involved.
Maintained swimming by individuals of D. gigas larger than 75 cm mantle length has been estimated as~0.5 m/s by archival tagging methods (Gilly et al. 2006;Stewart et al. 2012), but no comparable data are available for small individuals. Although fast swimming may not be necessary for prey capture, maintained swimming would be relevant to foraging, ultimately determining the volume of water in which a squid could search for aggregations of prey.
Prey accessibility could also explain diet variability between large and small squid if large squid forage deeper or have more success foraging at depth than do small squid. Vertical  (Robison 1972). Deeper habitats (400-500 m) are also utilized during the daytime by P. hoylei, which migrates into waters shallower than 100 m during the night (Young 1978). Archival tagging studies on large D. gigas in the GOC have revealed diel vertical migrations between nighttime near-surface (<75 m) and daytime mesopelagic (>300 m) depths (Gilly et al. 2006Bazzino et al. 2010). The vertical migratory capacity of small D. gigas is not well documented, but even if we assume the maximum depth range to be similar for large and small squid, maximum sighting distance increases with increasing eye size (Thomas et al. 2017). Large squid, with larger eyes, may have more success locating prey in darker, deeper habitats than small squid.

Prey availability is sensitive to environmental variability
Environmental conditions explained most of the variance in our diet dataset. Features of the physical habitat (SST, MEI) were important at longer lags than estimates of phytoplankton biomass (CHL) in explaining diet variability (~12 and 6 months, respectively, Table S3). This may reflect the importance of SST for recruitment of common prey species in the GOC, a transitional region accessible to tropical and subtropical species (Páez-Osuna et al. 2016), as well as the importance of food availability for post-recruitment retention of prey within a given habitat. Present environmental conditions (SST.01 and CHL.01) were only slightly less important than past conditions (MEI.13,SST.11,and CHL.06) in explaining variance in our CART analysis, suggesting that both recruitment and retention determined by past conditions as well as habitat selection under present conditions are important drivers of prey variability.
Groups of squid with similar diets determined using CART analysis also experienced similar trajectories of environmental conditions connecting past conditions to the present, and the habitat preferences of their main prey matched these intervening environmental conditions (summarized in Fig. 6). Squid in Nodes 12 and 13 experienced strong La Niña conditions throughout the year prior to capture and consumed a high proportion of Engraulis (Fig. 6B), a prey that prefers relatively cool, productive conditions associated with La Niña in the GOC (Kahru et al. 2004;Lluch-Cota et al. 2007;Checkley et al. 2009). B. panamense, a tropical, shelf-associated species (Robison 1972;Wisner 1974), was the dominant prey item in shelf habitats off Baja and Sonora when neutral ENSO conditions (MEI = −0.5 to 0.5) predominated throughout the year prior to capture (Fig. 6A) and during extended ENSO-positive, warm conditions that were relatively productive (Fig. 6C, D). When the GOC was warm and unproductive, large D. gigas fed on a high diversity of prey dominated by cephalopods regardless of intervening conditions (Fig. 6E). Under these same conditions, small squid fed mostly on crustaceans regardless of habitat (Fig. 6F).
Primary production in the GOC typically decreases during and immediately following El Niño events (Kahru et al. 2004), Fig. 6. Spatial trends of dominant prey for Dosidicus gigas in the GOC under all sampled environmental scenarios. Past MEI (MEI.13; "+" (positive) or "−" (negative)), ENSO phase during intervening 12 months ("Inter. ENSO") between past and present conditions ("La Niña" [MEI < −1.0], "neutral") [−1 < MEI < 1], "El Niño" [MEI > 1.0], or "ALL"), and present SST (SST.01; "warm," "cool," or "variable") are indicated above each map separated by vertical lines. Panels A and B illustrate dominant prey under neutral or relatively cool conditions. Panels C-F illustrate dominant prey under relatively warm conditions with variable chlorophyll a concentrations ("High CHL", or "Low CHL"). Past and present conditions were taken from the CART analysis and average intervening conditions were estimated from raw MEI values. The size classes of squid sampled under each set of conditions are represented as either large or small silhouettes. The nodes described are given in the bottom right corner of each panel. Blue indicates fish prey (dark blue-Benthosema, light blue-Engraulis), dark red indicates cephalopod prey, light green indicates crustacean prey, and white crosshatching indicates high prey diversity (D > 0.7).
which are known to induce changes in the oceanic environment of GOC, although local effects vary greatly among regions of the Gulf. Primary production in the GOC exhibits higher seasonal variability and relatively strong correlation with ENSO along the east coast (a wind-driven upwelling corridor) and southern Midriff Islands (a tidal mixing region), whereas lower variability and weaker correlation with ENSO occur along the western margin (Lluch-Cota et al. 2010). The impacts of this variability on middle trophic levels such as myctophids or squid are unknown. Similarly, decadal-tomultidecadal variability has seldom been investigated in the GOC due to the lack of long-term synoptic environmental data series (Lluch-Cota et al. 2010). However, since the 2009-2010 El Niño, wind-driven productivity in the central GOC has been anomalously low (Robinson et al. 2016). From 2014 to 2017, strongly ENSO-positive conditions prevailed in the GOC, and the 2015-2016 El Niño was followed by ENSO neutral conditions rather than a La Niña (Fig. S6). Squid in low-productivity conditions, regardless of size or capture region, fed on a relatively high diversity of prey and incorporated higher proportions of lower trophic-level prey (crustaceans, pteropods, and small cephalopods) than during highproductivity conditions (fishes; Ruiz-Cooley et al. 2006;Espinoza et al. 2017).
Biological and behavioral responses of prey to warming inferred from size-structure In addition to surface conditions and primary productivity affecting prey recruitment, prey availability can be modified by the effects of ocean temperature on body size and the vertical distribution of suitable habitats for different prey species (Cade and Benoit-Bird 2015;Cheung et al. 2013;Sheridan and Bickford 2011;Stramma et al. 2012). Many midwater species, regardless of vertical migratory behavior, exhibit ontogenetic descent, where larger individuals tend to occupy deeper daytime habitats than smaller individuals (Clarke 1973;Brodeur and Yamamura 2005). In general, mean prey size increased with increasing D. gigas mantle length for the three prey species measured, but the range of sizes of consumed prey remained similar across all squid sizes sampled. A consistent size range of consumed prey suggests that large and small squid are equally able to consume prey across the size range of prey sampled and that changes in mean prey size for these three small prey species reflect changes in the size-structure of prey populations.
B. panamense is a tropical species whose average size in squid stomachs did not appear to be dramatically affected by ENSO conditions. However, small squid captured during El Niño conditions exhibited the strongest positive relationship with B. panamense length, suggesting that large B. panamense were more accessible to these small squid than to squid captured under different environmental conditions. We also observed a trend of increasing P. hoylei size with increasing D. gigas size during El Niño conditions in regions with deeper bathymetry (Baja shelf and central pelagic). Compression of the vertical habitat used by B. panamense and P. hoylei due to surface warming during El Niño conditions and/or subsurface warming observed in the GOC since 2014 (Frawley et al. 2019) may have condensed groups of prey or reduced the distance between patches of prey, making larger prey more accessible. The negative relationship between B. panamense and D. gigas size during La Niña conditions could be explained by a change in diet to Engraulis, a shallower dwelling species than B. panamense. The average length of T. mexicanus decreased with increasing MEI and responded at a much longer lag than for either B. panamense or P. hoylei, which may be explained by its longer life span and utilization of deeper habitats (Robison 1972;Moser et al. 1984;Arkhipkin 2004). Abnormally warm surface temperatures would negatively affect the growth rate of juvenile T. mexicanus and result in a reduction in the average size of the population at depth following the ontogenetic descent (Jonsson and Jonsson 2014).
Specialization on relatively small, abundant prey D. gigas is considered a generalist predator and an ecological opportunist, taking advantage of new habitats and resources as they become available (Field et al. 2013). While it is clear that D. gigas is opportunistic, most notably during its invasions into waters of the Northern California Current (Field et al. 2007(Field et al. , 2013, this species exhibits similar or lower diet diversity than other oceanic predators of similar size (Jaquemet et al. 2011;Choy et al. 2013). Given the low variability in prey selectivity between individuals across a large size range and disparate habitats (Markaida and Sosa-Nishizaki 2003;Field et al. 2013;Alegre et al. 2014), we propose that D. gigas specializes on relatively small, abundant prey that form large aggregations. Squid, whose arms extend their functional gape, are less gape-limited than fish at a given body size and are able to capture relatively large prey (Rodhouse and Nigmatullin 1996); however, D. gigas larger than 35 cm mantle length (~1 kg) in the GOC consume smaller prey on average than similar sized fish predators (Scharf et al. 2000;Barnes et al. 2010). In fact, 10 kg D. gigas in this study still largely consumed 1-5 g myctophids.
A continued preference for small prey with increasing body size is likely greatly facilitated by the long, attenuate arm tips of D. gigas that have two to eight times more suckers on their arms than confamilial species (Roper et al. 2010). Fine arm tips allow squid of large body size to handle small prey, and D. gigas is by far the largest ommastrephid and the only family member with such fine arm tips (Roper et al. 2010). This arm morphology and feeding strategy may be directly responsible for the success of D. gigas throughout its range, because competition with similar sized fishes would be expected to be reduced with increasing body size through the continued consumption of small, relatively low trophic-level prey.

Conclusions
Over the past two decades, D. gigas has exhibited extreme variability in adult body size in the GOC but little variability in its trophic ecology with respect to body size. We provide an updated diet description for small individuals of D. gigas from the GOC and the most thorough time-series analysis of D. gigas diet to date. D. gigas exhibits a shift in prey selectivity at approximately 25 cm mantle length to favor higher trophic-level prey, but diet variability across the size-range sampled was largely driven by ENSO conditions and ocean temperature. Between 25 and 85 cm mantle length, D. gigas largely focuses on small, abundant prey. Feeding on large numbers of relatively small, high calorie prey in both relatively cool (anchovies) and relatively warm, productive conditions (myctophids) is likely necessary to support growth to large body sizes before maturation. When warm, unproductive conditions prevail, only small squid are present in the GOC and have diets dominated by euphausiids and pteropods, prey with relatively low caloric value. Shorter lived, shallower prey species responded to environmental variability at shorter temporal lags than longer lived, deeper prey species, suggesting that deeper dwelling species may be somewhat resilient to short-term variability in surface oceanographic conditions.
Both the life history of D. gigas and the composition of its forage respond to environmental variability in the GOC, but the relative contributions of temperature and food availability on the size of mature individuals remain unclear. From a prey availability standpoint, large squid that support productive fisheries can flourish in multiple environmental conditions by switching between the dominant small, energetically rich prey in different conditions. This variability may increase both the resilience of fisheries targeting D. gigas as well as the impacts on fisheries with which it overlaps. The observation that the environment was a stronger driver of variability in D. gigas trophic ecology than body size highlights the complexity of ecological interactions in a changing ocean. More effort to quantify the trophic response of other predators to environmental variability will certainly improve our understanding of and ability to effectively manage dynamic pelagic ecosystems.