Warmer, earlier, faster: Cumulative effects of Gulf of Alaska heatwaves on the early life history of Pacific cod

Anthropogenic climate change is pushing environmental conditions outside of the realm of historical variability (Williams and Jackson, 2007; Litzow et al., 2020; Litzow et al., 2021) at rates faster than many species can adapt (Sih et al., 2011; Crozier and Hutchings, 2014). The increase in frequency and intensity of marine heatwaves (MHWs) is an example of the rapid changes in environmental conditions caused by human-induced climate change (Laufkötter et al., 2020). Responses to MHWs have been tracked across trophic levels, and ecosystem, community, and population recovery after temperatures cool can last for years (Caputi et al., 2019; Arimitsu et al., 2021; Suryan et al., 2021). Yet our ability to predict future population responses requires an understanding of individual responses to temperature variation (Litzow et al., 2021). These individual responses include phenological and growth responses during extreme warming that can alter demographics and survival trajectories (Crozier and Hutchings, 2014; Reusch, 2014; Audzijonyte et al., 2016; Fiksen and Reglero, 2022). However, individual-level measurements spanning MHWs are difficult to quantify and rarely available.

Warmer temperatures can cause ectotherms to grow faster as juveniles, mature earlier, and reach a smaller maximum body size (the temperature-size rule [TSR]; Atkinson, 1994); however, responses to warming also differ across individuals, populations, and species. Species from higher trophic levels respond more variably to warming (Deutsch et al., 2022), and growth trajectories are more difficult to predict as conditions warm and food availability (prey quantity and quality) changes (Jobling, 1997; Jutfelt et al., 2021). Additionally, metabolism, growth, and maturation may differ depending on the number of generations exposed to warmer temperatures (Wootton et al., 2022). The lack of a universal mechanistic explanation of the TSR makes unclear if processes related to metabolism or to life history optimization drive TSR patterns (Audzijonyte et al., 2022). Without more definitive mechanistic understanding, phenological and growth responses to the extreme but relatively brief warming associated with MHWs remain difficult to anticipate.

Predicting responses to warming also requires disentangling the interwoven relationships between phenology, age, growth, and size, which requires detailed information on individuals. Empirical determination of age and growth are needed to separate ontogenetic and environmental effects and to quantify the effects of age versus growth on size (Morrongiello et al., 2014; Morrongiello and Thresher, 2015; Jørgensen et al., 2020). Additionally, individual responses to temperature can be influenced by biological and oceanographic variability (Koussoroplis et al., 2014; Dahlke et al., 2017; Thalmann et al., 2020), prior experience (Hurst et al., 2012; Herfindal et al., 2015), and genetics (Meffe et al., 1995; Burt et al., 2011). Quantifying the relationships among phenology, age, growth, and size is the first step for understanding the effects of temperature on individuals and cumulative impacts on populations.

The effect of temperature on phenology, growth, and size may be particularly important in early life stages. Early life stages are typically more sensitive to changes in environmental conditions than mature life stages (Rijnsdorp et al., 2009), as their smaller size and earlier developmental stage results in lower total energetic reserves (Post and Parkinson, 2001), greater susceptibility to predation (Miller et al., 1988; Paradis et al., 1996), and smaller geographic ranges (Ciannelli et al., 2022). Early life experiences can also establish legacies, such as carryover effects or compensatory growth, that influence phenotypes through later life stages and affect how individuals respond to later environmental conditions (Hector and Nakagawa, 2012; Almeida et al., 2021). Thus, the responses of individuals in early life stages can provide advance warning of a cohort’s recruitment and long-term performance (Morrongiello et al., 2014; Laurel et al., 2017; Litzow et al., 2022; Sakamoto et al., 2022).

Recent MHWs in the Gulf of Alaska (2014–2016 and 2019) dramatically altered the ecosystem and directly affected all life stages of many species due to their prolonged duration or intensity (Arimitsu et al., 2021; Suryan et al., 2021). These MHWs had dramatic negative consequences for Pacific cod (Gadus macrocephalus), which experienced a steep population decline and fishery closure (Barbeaux et al., 2020a) potentially due to an increase in metabolic demand and decrease in prey (Barbeaux et al., 2020b). The Gulf of Alaska MHWs were also linked to reductions in spawning habitat (Laurel and Rogers, 2020), temporal overlap between cod larvae and their prey (Laurel et al., 2021), and recruitment to coastal nursery habitats in the summer (Abookire et al., 2022). Despite these negative effects, pre-recruits were larger and in better condition in early summer (Abookire et al., 2022). However, age, growth, and size must be disentangled to determine if these size and condition responses were due to faster growth and/or earlier hatching (Jørgensen et al., 2020).

To quantify the relative importance of temperature on hatch phenology, size-at-age, and early growth of juvenile Pacific cod collected in coastal nurseries, we combined field collections, otolith structural analysis, and several statistical approaches. We also determined if, and how, these metrics have changed since the recent MHWs. We expected earlier hatch dates and larger sizes-at-age since the start of MHWs at least partially due to faster incubation and growth rates due to warmer temperatures, but we anticipated that the underlying relationships between temperature, phenology, and growth would be similar in years before, during, and between MHWs. We discuss mechanisms driving individual responses to temperature, which will be crucial to understand population and community responses to warming and to identify climate-ready fisheries.

Pacific cod inhabit all parts of the water column during their early life history (Laurel et al., 2023). They spawn in deep waters (75–200 m; Ketchen, 1961; Stark, 2007; Neidetcher et al., 2014) during January and February (Laurel et al., 2023). Eggs hatch near the sea floor (Bian et al., 2014) before quickly moving toward surface waters in the spring (March–April; Hurst et al., 2009). Pelagic larvae feed on zooplankton for 2–3 months (Laurel et al., 2023), with pelagic juveniles (20–40 mm standard length [SL]) observed as late as June (Dunn and Matarese, 1987). Mid-July is considered the settlement period for gadids near Kodiak because that is the earliest that they are captured in nearshore seines (BA Knoth and BJ Laurel, unpublished data).

The Kodiak Island Beach Seine Survey conducted by the National Oceanic and Atmospheric Administration (NOAA) Alaska Fisheries Science Center (AFSC) has collected age-0 Pacific cod during annual July surveys since 2006 in nearshore nursery areas near Kodiak Island (Anton Larsen Bay and Cook Bay; Figure 1). We used available archived fish stored at –20°C and collected during 2007, 2009, 2010, and 2012–2019 in our analyses (Tables 1 and S1). Sampling occurred at 16 sites spread evenly across Anton Larsen Bay and Cook Bay within a 1-week period beginning mid-July (range: July 6–25). Surface salinity and temperature (<2°C) were similar in the two bays on sampling days, and annual catch per unit efforts were similar to previous studies (Laurel et al., 2016). At each site, a 36 m demersal beach seine was deployed from a boat and pulled to shore by two people standing a fixed distance apart. The seine was 1 m deep with 13 mm mesh at the wings, expanding to 2.25 m in the middle with a 5 mm delta mesh in a bag-end. Seine sites were 2–4 m deep and were mostly sampled within 2 h of low tide (Laurel et al., 2007). Juvenile Pacific cod within each haul were identified, counted, measured (total length, mm), and then frozen at –20°C. As an indicator of coastwide annual abundance, we used estimates of catch per unit effort (no. fish set⁻¹) from a zero-inflated negative binomial Bayesian regression model that included the date of sampling as a smoothed non-parametric term and the nested site within bay as a group-level random effect (Table 1; Litzow et al., 2022).

We subsampled the juvenile collections for otolith analysis. We either examined all fish or, when more than 20 juveniles were collected, we randomly selected individuals from all vegetated sites (eelgrass and macroalgae), which were where fish were most consistently collected across years. Both sagittae (the largest of the three otolith pairs) were extracted from defrosted juvenile Pacific cod that had been weighed (to 0.01 mg) and measured (SL, to 0.05 mm). We prepared transverse sections of the left or right sagittae by mounting the otoliths on a glass microscope slide with thermoplastic resin, polishing the otolith on the posterior side, flipping the otolith, and continuing to polish until we exposed the core region. We used 3M^TM Wetordry^TM paper (800–2,000 grit), Buehler® lapping film (3–30 µm grit), and alumina slurry (0.3 µm) during polishing. Polished sections were imaged with a compound microscope at 400× magnification. Hatch check size and daily periodicity in juvenile Pacific cod otoliths have been previously confirmed (Narimatsu et al., 2007); therefore, we counted and measured the width of daily increments along the proximal–distal axis to estimate age and growth for each individual. Each otolith was interpreted by 1–3 readers 1–2 times using ImagePro Premier®, ensuring that each otolith was read independently 2–3 times. If independent counts varied by >10%, otoliths were revisited and discrepancies resolved. The read used to provide ages and increment widths in the analyses was randomly selected from multiple reads with <10% variation.

To develop a thermal history for each fish and to identify MHWs, we interpolated monthly temperature data from the oceanographic station GAK1 (http://research.cfos.uaf.edu/gak1/data/TimeSeries/). Juvenile Pacific cod captured at Kodiak Island are thought to have hatched within a geographic area represented by the GAK1 data and to settle in nurseries close to where they are hatched (Cunningham et al., 2009; Hinckley et al., 2019). GAK1 provides the longest and most complete seasonal temperature record for the Gulf of Alaska and reflects broader regional temperatures (see figure S2 in Laurel and Rogers, 2020). We further confirmed the applicability of the GAK1 data to temperatures experienced by early life stages of Pacific cod in the region by comparing GAK1 temperatures to those collected in Trident Basin, AK, over our study period, which provides a record of thermal variability in the Kodiak Island nearshore (Figure S1). These temperatures were significantly and highly correlated (ρ = 0.95, p < 0.001; Figure S2), similarly to what was documented in Laurel and Rogers (2020) with temperatures from the Shelik of Strait during the early larval period. Therefore, we assumed that GAK1 temperatures were representative of the likely experiences of Pacific cod eggs, larvae, and juveniles within the region.

Given that Pacific cod early life history involves hatching in demersal waters before moving up the water column (see above), we interpolated the vertical GAK1 monthly measurements to derive daily temperatures at depths of 100–250 m (hereafter “bottom”) as well as depths of 0–30 m (hereafter “surface”) using a cubic spline. We assigned a mean temperature value to each day of a fish’s life. For 0–5 days post-hatch (dph), we assumed individuals experienced bottom temperatures. After 5 dph, we assumed individuals experienced surface temperatures because they achieve neutral buoyancy around 3 dph, have strong surface orientation, and could feasibly swim to the surface within 1–2 days (Hurst et al., 2009; Li et al., 2015).

MHWs were identified as anomalous warm events during which temperatures were warmer than the 90th percentile of the 30-year baseline period of 1983–2012 for 5 days or more (Hobday et al., 2016; Hobday et al., 2018). We also interpolated GAK1 temperature values for this analysis. All heatwaves were identified and categorized using the functions “ts2clm,” “detect_events,” and “category” in R package “heatwaveR” (Schlegel and Smit, 2018), specifying the climatological period of January 1, 1983, to December 31, 2012, and that the data were collected in the northern hemisphere.

Although our analyses focused on evaluating the relationships between temperature and Pacific cod phenology, size, and growth, we binned the years into three categories, before (2007, 2009, 2010, 2012–2014), during (2015–2016, 2019), and between MHWs (2017, 2018), to determine if temperature relationships differed between these three periods. While heatwaves did occur in surface waters before 2015, all but one were moderate heatwaves (the mildest category) and occurred in July or later, which was outside of or toward the end of the growth period in our analyses. No heatwaves occurred before July until 2015 (Table S2, Figure S3). The prolonged heatwave colloquially termed “The Blob” (Bond et al., 2015) began in 2014 (November 27, 2014, to April 15, 2015; Table S2); however, we did not categorize 2014 as a MHW year because, based on the GAK temperature data, the heatwave started after the 2014 fish had been collected. In other words, it overlapped with the spawning, hatching, early life, and growth periods of juveniles starting in 2015, not in 2014. We included 2017 and 2018 in their own category, as these years were warmer than the before MHW period (Figure 2; Suryan et al., 2021). MHWs reoccurred during 2019, demonstrating that more frequent, severe MHWs may become typical (Laufkötter et al., 2020).

We estimated hatch date distributions for the full cohort by determining daily age from sampled otoliths (n = 8–33 per year) and applying annual age–length relationships to the remaining catch (see Table S3 and Figures S4 and S5 in Text S1). To determine if Pacific cod hatch dates differed among MHW periods, we used a linear mixed model (R package “lme4”; Bates et al., 2015) with the response variable of individual hatch date, fixed categorical effect of heatwave designation (before, during, or between MHW), and a random intercept of year to incorporate the dependency among individuals within the same year and to account for any unexplained variation in annual conditions. We followed the linear mixed model with a type-II Wald χ² test (R package “car”; Fox and Weisberg, 2019) to evaluate significance and computed the estimated marginal mean (least-squares means) hatch dates for each period (R package “emmeans”; Lenth, 2021).

We quantified the direct effect of incubation temperature on hatch phenology. Pacific cod embryonic development requires approximately 110.5 cumulative degree days (CDD; i.e., the sum of mean daily temperatures in °C; Narimatsu et al., 2007; Laurel et al., 2008). For each fish, we determined the number of days before their hatch date required to reach 110.5 CDD and averaged the bottom temperatures during this incubation period. We then determined how much variation in hatch dates could be attributed to incubation temperature. To make this determination, we compared the observed shifts in hatch dates across years (using the overall marginal means from the linear mixed model described above) with the estimated difference in incubation duration based on bottom temperature (proportion difference in hatch dates due to incubation temperature = difference in calculated incubation duration/observed difference in mean hatch dates). This approach allowed us to determine if changes in incubation duration due to warming could explain the observed shifts in hatch dates.

For all statistical analyses herein, we validated the models by examining plots of residuals versus fitted values, residuals versus predictor variables, and quantile–quantile plots. If these plots indicated that model assumptions were not met, we modified the family of analysis, log-transformed the response variable, or determined if a higher polynomial was more appropriate for one of the predictor variables and specified this change in the relevant text. All analyses were conducted in R (R Core Team, 2021), and we visualized our results throughout the present study with R package “ggplot2” (Wickham, 2016).

Size-at-age was estimated as SL (mm) of fish at capture with ages either obtained from otoliths or applied from year-specific age-length relationships. To determine if Pacific cod size-at-age differed for individuals hatched before, during, and between MHWs, we used a linear mixed model (R package “lme4”; Bates et al., 2015) with the response variable of individual SL, the fixed main and interactive effects of heatwave designation (categorized as before, during, or between MHW) and age (continuous, scaled and centered), and a random intercept of year. We followed the linear mixed model with a type-III Wald χ² test (R package “car”; Fox and Weisberg, 2019) to evaluate significance. To aid in visualization of the results, we extracted the predicted main and interactive effects, including predicted size-at-age for periods before, during, and between MHW (R package “effects”; Fox and Weisberg, 2019).

To better understand the mechanisms driving differences in size-at-age, we evaluated daily growth responses to temperature, hatch date, MHW designation, and age. Daily growth rates were estimated with back-calculated size-at-age (mm) from otolith increment widths. The relationship between SL and otolith radius was close enough to linear (see Table S4 and Figures S6 and S7 in Text S2) that results from the modified Fry and linear biological intercept were similar (Vigliola et al., 2000; Vigliola and Meekan, 2009). Therefore, we used the biological intercept back-calculation model (Campana, 1990) with published values for Pacific cod biological intercepts (Narimatsu et al., 2007). The length of an individual at age a (L_a) is:

where L_c is the length at capture, O_a is the otolith radius at age a, O_c is the otolith radius at capture, and L₀ and O₀ define the biological intercept of the length (L) and otolith radius (O) at hatching (Figure S7).

We assessed the relationships between daily growth and temperature as well as MHW designation and age using generalized additive mixed models (GAMM, R package “mgcv”; Wood, 2017). The GAMM allowed us to avoid preemptively assume the shape of the relationship between predictor variables and growth and to include random effects to incorporate the dependency among daily growth observations from the same fish and year, as well as to account for ontogeny. We modeled the response variable of back-calculated daily growth rates (mm d⁻¹) as a function of MHW designation, age, and temperature:

We natural-log-transformed daily growth (G) to improve model diagnostics. Our predictors included a main effect of the categorical variable MHW (HW, levels before, during, and between) and the three-dimensional tensor product smooth function (te)of temperature (T), daily age (A), and MHW (HW). The tensor product smooth-evaluated the effect of temperature at different ages on growth while allowing these relationships to vary by MHW designation. All smooths were implemented with cubic regression splines. Fish ID (unique identifier for each fish) and Year were random intercepts, and A was a random slope for the random intercept ID to account for ontogenetic effects (Morrongiello and Thresher, 2015). We used restricted maximum likelihood and scaled and centered all continuous predictor variables (A, T). We confirmed that concurvity among variables was reasonable for the model (see Tables S5–S8 in Text S3). To determine how size affected growth responses, we ran an identical analysis with relative daily growth (mm mm⁻¹ d⁻¹), which we also natural-log-transformed to improve model diagnostics (see Table S9 and Figure S8 in Text S4).

Our heatwave designations generally grouped cooler years as “before MHWs” and warmer years as “during MHWs,” but inter- and intra-annual variation in bottom and surface temperatures demonstrated that a broad range of temperatures was experienced in all three periods (Figure 2A–C). Incubation and bottom temperatures were consistently warmer since MHWs, although 2017 and 2018 (“between MHWs”) had temperatures similar to the warmest years before MHWs (Figure 2D and E). Surface temperatures demonstrated greater variation within years, ranging from 2.5°C to >12.0°C between late March and August, than bottom temperatures. The warmest annual mean surface temperature occurred during 2015 and coolest during 2012 (Figure 2F).

Compared to before MHWs, hatch dates occurred 14 days earlier on average during MHWs and 26 days earlier between MHWs in the Gulf of Alaska (Tables 2 and 3; Figure 3). The latest annual mean hatch dates occurred during the before years 2009 and 2013 (April 17 in both years), which was 43 days later than the earliest annual mean hatch date (March 5) in 2018 (Table S1). The difference in average incubation temperatures before and during MHWs was approximately 1.8°C, which corresponded with an approximate 7.5-day difference in hatch dates. Therefore, of the estimated 14-day shift in hatch dates during MHWs, approximately 53% could be attributed to faster development due to increases in incubation temperature. The difference in average incubation temperatures before and between MHWs was approximately 1.0°C, which indicated that only about 4 days of the 26-day difference in hatch dates, or 16%, could be attributed to faster development due to increases in incubation temperature. Although incubation temperatures between MHWs were an average of 0.8°C cooler than those during MHWs, hatch dates were approximately 12 days earlier between MHWs.

Size-at-age differed among MHW periods. There was a reversal in the pattern of size-at-age between MHW status that resulted in an interaction between age and MHW (Table 2). Early in life, fish were similar in size-at-age across MHW periods, with a trend of lower average size-at-age between MHWs (Figure 4B). However, after approximately 120 days, fish were larger-at-age during and between MHWs than before MHWs. This pattern resulted in a mean difference in size of 11.4 and 7.1 mm at 132 days between the periods before and between MHWs and the periods before and during MHWs, respectively (linear mixed model results; Tables 2 and 3; Figure 4B). Size-at-age during MHWs was generally larger than before MHWs for all ages, with model results not accounting for age indicating that individuals during MHWs were approximately 4 mm larger than before or between MHWs.

Growth-at-age responded to temperature differently among the MHW periods (GAMM results: adjusted R² = 0.71; Table 4, Figures 5 and 6). Early (<14 dph) growth rates were generally slower during MHWs, with slower growth at all temperatures experienced at 10 dph (Figure 6), although higher temperatures were experienced during MHWs than in the other periods. From 25 to 50 dph, growth increased faster with age during MHWs (note narrowness of model-predicted contours in Figure 5B). Interestingly, growth was slowest at intermediate temperatures during all three periods at 40 dph, with the temperature considered “intermediate” increasing by approximately1°C from before to between to during MHWs (GAMM results: Figures 5 and 6E, F). By about 60 dph, growth rates increased primarily due to temperatures rather than age in all MHW periods, with higher growth rates achieved at higher temperatures during and between MHWs than before (GAMM results: Figures 5 and 6G). Model-predicted growth during and between MHWs out-paced growth before MHWs above about 7.5°C by approximately 50 dph (Figure 5). Some of the absolute growth differences before and between MHW periods after 60 dph were due to larger body size between MHWs given similar relative growth rates during that period (note greater overlap of 95% confidence intervals of modeled data in Figure 6H compared to Figure 6G; Figure S8). However, accounting for body size increased the range of temperatures for which there were differences in growth during MHWs where relative growth was slower below approximately 8°C but higher above 10°C compared to before and between MHWs (Figure 6H).

Our results demonstrate how Pacific cod juvenile hatch phenology, size-at-age, and growth differed before, during, and between the periods of recent prolonged warming in the Gulf of Alaska. For juveniles that survived to settle in coastal nurseries, we observed much earlier (approximately 2–4 weeks) hatch phenology and slightly faster growth since 2015, with patterns persisting between MHWs despite a return of cooler temperatures. Individuals were also larger at age during and between MHWs, but only for older individuals (>120 days old); observed increases in growth alone could not account for those larger sizes. Temperature was related to the observed differences in phenology and growth, often in the ways predicted by metabolic theory such as accelerated egg incubation rates leading to earlier hatch dates (Bruno et al., 2015). Approximately 53% of the shift to earlier hatch dates during MHWs could be explained parsimoniously by temperature-dependent egg incubation rates, indicating that about 47% of the change in phenology during MHWs was due to indirect effects of temperature such as parental behavior (Schneider et al., 2010; Cannaday and Farmer, 2022), changes in larval source or selective survival (Hinckley et al., 2019; Fiksen and Reglero, 2022), and/or co-occurring environmental changes. However, between MHWs, much more of the change in phenology was due to other factors, as only 16% of the shift in hatch dates could be explained by faster incubation rates.

The relationships between temperature and growth were not straightforward, most obviously with the change in growth-at-age with temperature among the MHW periods. Similar temperatures could produce differences in growth rates >0.1 mm d⁻¹ at the same age among periods. While some of these differences were related to size, differences in size-specific growth at the same age and temperature were also significant. Thus, temperature alone cannot explain the observed growth differences between MHW periods. Individual growth responses to temperature before, during, and between MHWs were likely also affected by environmental changes associated with MHWs, such as prey abundance, timing, energetic quality, and composition (Arimitsu et al., 2021; Suryan et al., 2021); previous experiences (Hurst et al., 2012; Fuiman and Perez, 2015); and/or population composition, such as shifts in spawning populations (Rogers and Dougherty, 2019) or selection (Murphy et al., 2014). By quantifying the relationships between temperature, phenology, age, and growth, we provide valuable knowledge about the interrelated way in which warming temperatures alter early life history and could contribute to larger population consequences through survival or carryover effects.

Our detailed examination of wild fish phenology and daily growth spanning a broad range of thermal conditions also contributes to understanding broader-scale responses to warming. Mechanistic explanations for population-level responses to warming, such as the TSR, predict faster juvenile growth rates and larger sizes (Audzijonyte et al., 2022; Deutsch et al., 2022; Wootton et al., 2022); however, these explanations do not fully encapsulate our key observations. Importantly, neither metabolic nor life history optimization mechanisms for the TSR predict the similar size-at-age and faster growth-at-age during early life (<14 dph) in the cooler before MHW period. Additionally, the eventual approximately 11 mm larger size-at-age of juveniles at 132 dph between MHWs cannot be wholly explained by differences in growth rates. Daily growth during and between MHWs surpassed growth before MHWs at all temperatures by 40 dph, with a mean difference in growth of 0.04–0.05 mm d⁻¹ after that rate. These faster growth rates only account for about 3.7–4.7 mm between 40 dph and 132 dph, which does not explain the size differences even if we were to assume equivalent starting sizes at 40 dph. Therefore, positive size selection could have contributed to these differences in size-at-age (Murphy et al., 2014).

Proposed TSR mechanisms do not explicitly consider phenological shifts associated with warming that create age variation within life stages. More of the increase in juvenile size-at-capture since 2015 can be explained by earlier hatch dates and older ages than by increased growth rates. Hatch dates were 14 days and 26 days earlier on average during and between MHWs, respectively, and up to 43 days earlier in pairwise comparisons of years. Additionally, observed increases in size-at-age since 2015 cannot be explained by the growth rates observed, highlighting the likely role of enhanced size-selective mortality during warmer times. Mean growth rates after 75 dph of 0.44 mm d⁻¹and 0.46 mm d⁻¹ could lead to fish that were an average of 6.2 mm and 12.0 mm longer, respectively, by July or up to 19.8 mm longer in 2018 only due to their older ages. Such a large size difference could affect competition (Van Allen and Rudolf, 2016), over-winter survival (Post and Evans, 1989; Litz et al., 2017), and “carry over” to older ages (Chelgren et al., 2006; Almeida et al., 2023). The relative importance of hatch phenology and growth to our results highlights the need to consider both, as well as size-selective mortality, in regard to life-history optimization under different thermal conditions and raises questions about how age variation fits into the mechanistic explanations of the TSR.

Although earlier hatching is a documented response to climate change (Edwards and Richardson, 2004; Menzel et al., 2006; Poloczanska et al., 2016), a shift of the magnitude we observed since MHWs is large. Earlier hatching may result in competitive advantages for larger, older individuals and/or an increase in predation on smaller, younger competitors (Rasmussen et al., 2014). Intra-specifically, these competitive advantages could further exaggerate the phenological shift of the survivors of the annual cohort. However, a phenological shift of the magnitude we observed likely caused Pacific cod juveniles to experience different early-life environmental conditions in each MHW period, which may have resulted in high mortality and contributed to the low July abundances during MHW years (2015, 2016, 2019).

A contributing factor to the shift in hatch dates beyond the direct effects of temperature could be selective survival due to changes in food, transport, and/or predation. Match–mismatch dynamics of larvae and prey can cause hatch date selection (Cushing, 1969, 1990; Durant et al., 2007; Peck et al., 2012; Neuheimer et al., 2018), with similarly large shifts in optimal hatch dates predicted from models for Atlantic bluefin tuna (Thunnus thynnus) during a MHW (Fiksen and Reglero, 2022). Mismatches and resulting mortality may be more likely to occur during MHW years because warmer temperatures reduce the period of time that larvae can rely on yolk sac reserves prior to first feeding (Laurel et al., 2008; Laurel et al., 2011; Laurel et al., 2021; Fiksen and Reglero, 2022). However, prey phenology may also shift with hatch phenology, reducing the period of mismatch. In the Gulf of Alaska, primary production timing shifted approximately 20 days earlier since MHWs (Laurel et al., 2021). More detailed data on primary and secondary production could be used to evaluate potential mismatches and may provide insight on the degree to which selection affected juvenile responses to MHWs.

Changes in behavior, demographics, or the source of spawners also could have contributed to the shift in hatch dates during and between MHWs. Warmer winters and springs cause earlier spawning in many fish species (Kjesbu et al., 2010; Schneider et al., 2010; Rogers and Dougherty, 2019; Cannaday and Farmer, 2022). Additionally, older and larger fish can spawn earlier (Kjesbu et al., 2010; Rogers and Dougherty, 2019) or over a longer period (Ohlberger et al., 2014) and provision their eggs to a greater extent (Berkeley et al., 2004). Thus, changes in the age or size structure of the adult population due to variable survival during MHWs could lead to changes in hatch dates during and between MHWs (Barbeaux et al., 2020b). Additionally, the population breadth or source of spawners may have changed since MHWs. The most likely spawning areas for juveniles that ultimately recruit to nursery grounds in northern Kodiak Island are 20–200 m depths (Dunn and Matarese, 1987; Cunningham et al., 2009; Hinckley et al., 2019); however, after MHWs, only individuals that had been spawned at shallower depths may have survived to the juvenile life stage due to thermal restrictions on spawning habitat suitability and the potential for cooler shallower waters (Laurel and Rogers, 2020). Given the potential for parental genetics and experience to influence offspring responses to temperature beyond hatching (Burt et al., 2011), determining if parental behavior, demographics, or composition changed among MHW periods is an important step in understanding population-level consequences of MHWs.

Size-at-age and growth rates differed among MHW periods, with periods of faster growth-at-temperature during and, especially, between MHWs, providing a carryover effect of larger sizes-at-age. Our daily growth results also indicated that surviving individuals responded to temperature differently among MHW periods. Although fish may not have experienced the temperatures we assigned given the range of potential spawning habitat (Hinckley et al., 2019), our temperatures are similar to other records in the region (Figures S1 and S2; Laurel and Rogers, 2020) and the relative change in growth due to changes in temperature should be reliable. Overall, while the differences in growth among MHW periods contributed to the differences in size-at-age, they could not account for the observed increases; therefore, identifying other factors regulating size will help determine the way in which MHWs affected early life stages of Pacific cod. Multiple mechanisms, such as food limitation, prior experience, and genetic (trait) selection, likely affected size-at-age and growth simultaneously.

Food limitation (quantity or quality) may explain variation in growth in relation to temperature and age. During the first few weeks of life, the lower growth rates and shallower increase in growth as temperatures warmed during MHWs may indicate that these individuals were unable to maximize growth due to higher metabolic demands and inadequate prey quantity or quality (Jobling, 1997; Hurst et al., 2010; Laurel et al., 2011; Fiksen and Reglero, 2022; Sakamoto et al., 2022). The shifts in timing of phytoplankton blooms (Laurel et al., 2021) and potential mismatch with prey may have affected growth as well as phenology. General trends in Gulf of Alaska warm and cool water zooplankton abundance were positive or neutral during the 2015–2016 MHW years (Suryan et al., 2021), but the composition of spring zooplankton communities in parts of the Gulf of Alaska differed during MHWs, with greater dominance of smaller, less nutritious copepods (Arimitsu et al., 2021). While the spring zooplankton communities appear to differ in warmer versus cooler years, this difference may not account appropriately for changes in timing, size, and developmental stage. Further evaluations of fish diet and energetic data would provide stronger evidence for the effects of changing zooplankton communities on fish early life stages.

Although food quantity and quality may have changed since MHWs, the relationship between temperature and growth-at-age also differed among MHW periods potentially due to plastic or genetic differences in the individuals that composed annual cohorts. Plastic changes caused by previous experiences of the juveniles or their parents can affect how individuals respond to environmental conditions throughout life (Burt et al., 2011; Burton and Metcalfe, 2014; Almeida et al., 2021). For example, if earlier faster growth allows individuals to achieve larger sizes, individuals may shift from zooplanktivory to consuming macroinvertebrates or fish at an earlier age, which can have cascading consequences in growth and development (Vigliola and Meekan, 2002; Fedewa et al., 2017; Jørgensen et al., 2020). Maternal contributions to eggs and experiences during the embryonic and early larval stages can cause phenotypic plasticity in temperature-dependent growth or epigenetic effects due to metabolic programing or environmental conditioning (Zambonino-Infante et al., 2013; Fuiman and Perez, 2015; Jørgensen et al., 2020); previous studies have found support for plastic temperature-growth responses in Gulf of Alaska Pacific cod populations (Hurst et al., 2012). Genetic differences in annual cohort composition could occur due to selection for specific traits (Meffe et al., 1995; Croisetière et al., 2010), which is a key component to climate adaptation (Crozier and Hutchings, 2014; Reusch, 2014). The low juvenile abundances in 2015, 2016, and 2019 indicate that mortality was high and selection could have been particularly strong in these years. Changes in parental source could also affect the genetic composition of the annual cohort, which may have occurred due to variation in currents (Weingartner et al., 2009; Hinckley et al., 2019) or selection of spawners (Rogers and Dougherty, 2019). Otolith δ¹⁸O isotopes could be used to verify thermal exposure (Kastelle et al., 2022), with core values potentially indicating if hatching location differed between cohorts; values through the larval and juvenile period could determine if thermal experiences differed. Overall, the different growth-at-temperature responses among MHW periods were likely caused by a combination of other environmental conditions and previous experiences and reflect the surviving population composition. To evaluate the potential for these different types of legacy effects, a variety of laboratory experiments and genetic information would be required, but these efforts would be worthwhile to better understand recruitment variability (Laurel et al., 2023).

MHW designations will likely become more arbitrary as thermal baselines shift and species and ecosystems reorganize to novel environments. For example, our decision to separate those years that were warmer but did not meet MHW criteria (2017, 2018) into the between-MHW period was based on thermal trends but also demonstrated a longer lasting influence of MHW conditions on Pacific cod biology. During MHW years (2015, 2016, 2019), earlier hatch phenology and faster growth was related to temperature, but the low abundance of juveniles captured during those years (Table 1; Abookire et al., 2021) provides support for poor survival and potentially strong selection as well. In 2017 and 2018, which were cooler but still warm, the earlier hatch phenology and faster growth persisted, with 2018 having the earliest hatch dates of all years in our analysis and the fastest growth and largest size-at-age consistently in the between-MHW period. However, juvenile abundance was high (Table 1; Abookire et al., 2021), potentially indicating that factors other than selective mortality were contributing to the maintenance of these MHW patterns in non-MHW years. Changes associated with community and ecosystem responses in the Gulf of Alaska also lingered during 2017 and 2018 (Suryan et al., 2021). Overall, the persistent evidence of MHW influences on Pacific cod juveniles in 2017 and 2018 demonstrates the value of focusing on thermal trends rather than binary categories of MHWs but also suggests that the effects of the most recent MHWs in the Gulf of Alaska may have lastingly altered populations.

We disentangled the complex effects of temperature on phenology, growth, and size, demonstrating that the responses of juvenile Pacific cod to temperature differed before, during, and between MHWs. During and between MHWs, individuals hatched earlier and were larger-at-age when they reached coastal nurseries than before MHWs. The relationships between temperature and hatch timing, size-at-age, and growth changed, however, during MHWs and persisted despite cooler temperatures between MHWs. Thus, our results highlight the need to understand how factors other than temperature interact to influence Pacific cod during their early life stages. Other environmental conditions (such as the prey community), previous experiences (epigenetic effects or carryover effects), size- or growth-based selection, or population composition (genetics) likely also influenced the observed responses of Pacific cod to the MHWs. Although the Gulf of Alaska does not appear to have shifted to a new ecosystem state based on time-series analyses to 2019 (Litzow et al., 2020), the MHWs may have been impactful enough to change the way in which Pacific cod populations respond to thermal conditions. The persistence of these new relationships will be important to follow. By improving understanding of the direct thermal effects of MHWs on early life stages of Pacific cod, our findings will be useful for explaining and predicting population responses to climate change (Barbeaux et al., 2020b; Litzow et al., 2021; Peterson Williams et al., 2021; Litzow et al., 2022) and evaluating mechanisms driving TSR patterns. While considering whole ecosystem changes caused by unprecedented and rapid climate change can be untenable, building understanding of the mechanisms driving individual responses to thermal conditions can improve our ability to anticipate population changes and manage species.

Data and code are available in a GitHub repository: https://github.com/LZoe/warmer-earlier-faster.

The supplemental files for this article can be found as follows:

Table S1. Information from age-0 juveniles used in the analyses.

Figure S1. Temperature comparison between GAK1 and Trident Basin temperature station.

Figure S2. Correlation between GAK1 surface temperatures and Trident Basin temperatures.

Table S2. Identified marine heatwaves (MHWs) during the 11 years juveniles were captured.

Figure S3. Identified marine heatwaves (MHWs) since 2005 including the 11 years juveniles were captured.

Text S1. Age-length key details.

Table S3. Type III SS ANOVA results evaluating the effect of SL and year on age.

Figure S4. Relationship between standard length (SL) and age of juveniles from our aged subsample.

Figure S5. Distributions of hatch dates as indicated with density.

Text S2. Fish length and otolith radius relationship.

Table S4. Type III SS ANOVA results evaluating the effect of radius, year, life stage on SL.

Figure S6. Linear relationships between otolith radius and standard length.

Figure S7. Chronologies of back-calculated absolute (mm d^–1) and relative (mm mm^–1 d^–1) growth.

Text S3. Evaluating concurvity.

Table S5. Overall concurvity values for each smooth in the absolute growth (mm d^–1) GAMM.

Table S6. Overall concurvity values for each smooth in the relative growth (mm mm^–1 d^–1) GAMM.

Table S7. Pairwise concurvity values for each smooth in the absolute growth (mm d^–1) GAMM.

Table S8. Pairwise concurvity values for each smooth in the relative growth (mm mm^–1 d^–1) GAMM.

Text S4. Generalized additive mixed model (GAMM) with relative growth (mm mm^–1 d^–1).

Table S9. GAMM results for the response of relative growth to temperature and age before, during, and between MHWs.

Figure S8. GAMM-predicted relationships between age-0 Pacific cod back-calculated relative growth (mm mm^–1d^–1) and temperature influenced by age.

The authors thank Thomas Murphy for assisting with the preparation of otolith samples and members of the NOAA AFSC field crews who collected juveniles through the Kodiak Island Beach Seine Survey. Thanks to John Cannaday for assistance with otolith interpretation. They also thank Dr. Lauren Rogers, Dr. Michael Litzow, the Miller lab, and members of the NOAA Fisheries Behavior Ecology Program for valuable feedback during the development of the analyses. Additional thanks to Drs. Esther Goldstein and Thomas Hurst for providing feedback on drafts of this manuscript.

This work was supported by the North Pacific Research Board, project No. 1903 (to JAM).

The authors declare there are no competing interests.

• Contributed to conception and design: LZA, BJL, HLT, JAM.

• Contributed to acquisition of data: LZA, BJL, HLT, JAM.

• Contributed to analysis and interpretation of data: LZA, BJL, HLT, JAM.

• Drafted and/or revised the article: LZA, BJL, HLT, JAM.

• Approved the submitted version for publication: LZA, BJL, HLT, JAM.

How to cite this article: Almeida, LZ, Laurel, BJ, Thalmann, HL, Miller, JA. 2024. Warmer, earlier, faster: Cumulative effects of Gulf of Alaska heatwaves on the early life history of Pacific cod. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2023.00050

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Maxime Geoffroy, Fisheries and Marine Institute of Memorial University of Newfoundland, St. John's, Canada

Knowledge Domain: Ocean Science