Warming climates are creating unprecedented environmental conditions, such as more frequent and intense marine heatwaves (MHWs), that directly impact phenology and growth of fish and other marine organisms. Understanding individual phenological and growth responses to temperature is critical to predict species and population responses to climate change; however, doing so requires disentangling the effects of temperature on phenology, size, and growth in wild populations. We quantified the relationships between temperature and hatch timing, size-at-age, and early growth in a population of Pacific cod (Gadus macrocephalus) affected by recent MHWs in the Gulf of Alaska. Pacific cod juveniles were collected near Kodiak Island, Alaska, USA, across 11 years, categorized as before (2007, 2009–2010, 2012–2014), during (2015, 2016, 2019), and between (2017, 2018) multiple recent MHWs. We estimated age and growth with otolith structural analysis. Hatching occurred on average 14 days earlier during MHWs and 26 days earlier between than before MHWs. Approximately 53% and 16% of these respective shifts in timing were attributable directly to warmer temperatures during incubation. Size-at-age was similar across periods at younger ages (90 days), but approximately 7 mm and 11 mm larger than before MHWs at older ages (132 days) during and between MHWs, respectively. These differences in size-at-age were partially related to growth responses that differed among MHW periods. However, observed differences in growth rate could not account for the observed increases in size-at-age. We found that temperature alone could not explain the changes in growth and phenology; thus, factors such as parental effects, epigenetics, and selection likely contributed. Our results indicate that spawn timing, size, and growth relationships based on historical thermal responses should be questioned in population forecasting as the global climate continues to warm.

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 early life history

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).

Biological collections and otolith analysis

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−1) 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).

Figure 1.

Juvenile Pacific cod collection sites near Kodiak Island, Alaska, USA. Age-0 juvenile Pacific cod were collected from nearshore nursery habitat in Anton Larsen Bay and Cook Bay (light blue triangles in inset) near Kodiak Island, Alaska, USA.

Figure 1.

Juvenile Pacific cod collection sites near Kodiak Island, Alaska, USA. Age-0 juvenile Pacific cod were collected from nearshore nursery habitat in Anton Larsen Bay and Cook Bay (light blue triangles in inset) near Kodiak Island, Alaska, USA.

Close modal
Table 1.

Information from age-0 juvenile Pacific cod used in the analyses

DesignationaYearN CapturedbN ProcessedcNo. Fish Set−1 dSL (mm)eAge (days)f
Before 2007 31 28 7.4 ± 2.9g 43.9 ± 0.5 98.7 ± 1.5 
2009 67h 33 26.6 ± 15.3 38.7 ± 0.5 87.5 ± 0.9 
2010 104 23 9.2 ± 4.8 45.2 ± 0.5 97.3 ± 1.3 
2012 46 33 171 ± 63.5 40.0 ± 0.7 103 ± 1.1 
2013 41 19 8.1 ± 3.5 49.4 ± 0.8 98.1 ± 1.5 
2014 38h 20 7.7 ± 4.2 50.4 ± 1.1 115 ± 1.8 
During 2015 1.0 ± 1.0 47.7 ± 1.3 98.9 ± 2.9 
2016 49 19 1.7 ± 0.8 63.5 ± 0.9 120 ± 1.3 
Between 2017 32 24 74.7 ± 24.3 59.9 ± 2.0 114 ± 2.3 
2018 46 20 104 ± 24.4 64.6 ± 2.0 133 ± 2.1 
During 2019 31 20 1.8 ± 1.1 49.4 ± 1.2 111 ± 1.8 
DesignationaYearN CapturedbN ProcessedcNo. Fish Set−1 dSL (mm)eAge (days)f
Before 2007 31 28 7.4 ± 2.9g 43.9 ± 0.5 98.7 ± 1.5 
2009 67h 33 26.6 ± 15.3 38.7 ± 0.5 87.5 ± 0.9 
2010 104 23 9.2 ± 4.8 45.2 ± 0.5 97.3 ± 1.3 
2012 46 33 171 ± 63.5 40.0 ± 0.7 103 ± 1.1 
2013 41 19 8.1 ± 3.5 49.4 ± 0.8 98.1 ± 1.5 
2014 38h 20 7.7 ± 4.2 50.4 ± 1.1 115 ± 1.8 
During 2015 1.0 ± 1.0 47.7 ± 1.3 98.9 ± 2.9 
2016 49 19 1.7 ± 0.8 63.5 ± 0.9 120 ± 1.3 
Between 2017 32 24 74.7 ± 24.3 59.9 ± 2.0 114 ± 2.3 
2018 46 20 104 ± 24.4 64.6 ± 2.0 133 ± 2.1 
During 2019 31 20 1.8 ± 1.1 49.4 ± 1.2 111 ± 1.8 

aMarine heatwave (MHW) designation.

bSample sizes of fish captured.

cSample sizes of fish from which otoliths were extracted and processed.

dEstimated abundance from a zero-inflated negative binomial Bayesian regression model (Litzow et al., 2022).

eStandard length-at-capture.

fApplied from an age-length regression.

gAll data with error bounds represent means with standard error.

hTwo fish were removed, one from each of these years, because they were far (>10 mm) outside the range of sizes of aged fish and thus our age-length key could not be applied to them reliably.

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 3MTM WetordryTM 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.

Temperature and MHW designations

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).

Figure 2.

Gulf of Alaska daily interpolated temperatures (A–C) and annual means ± SD (D–F). Daily interpolated bottom (100–250 m) temperatures are depicted for periods reflective of (A) the 4 weeks before the first hatch date (thin lines) and first-to-last hatch date (thick lines) and (B) the day the first individual hatched to 5 days after the last individual hatched each year. Daily interpolated surface (0–30 m) temperatures are depicted for (C) the period from 5 days after the first individual hatched until the day individuals were captured. Mean annual temperatures ± SD are presented for (D) averaged individual incubation temperatures, (E) interpolated temperatures in panel B, and (F) interpolated temperatures in panel C. Colors indicate the year in (A)–(C) and how each year was grouped in (D)–(F): before (black), between (blue), or during (light blue) marine heatwaves (MHWs).

Figure 2.

Gulf of Alaska daily interpolated temperatures (A–C) and annual means ± SD (D–F). Daily interpolated bottom (100–250 m) temperatures are depicted for periods reflective of (A) the 4 weeks before the first hatch date (thin lines) and first-to-last hatch date (thick lines) and (B) the day the first individual hatched to 5 days after the last individual hatched each year. Daily interpolated surface (0–30 m) temperatures are depicted for (C) the period from 5 days after the first individual hatched until the day individuals were captured. Mean annual temperatures ± SD are presented for (D) averaged individual incubation temperatures, (E) interpolated temperatures in panel B, and (F) interpolated temperatures in panel C. Colors indicate the year in (A)–(C) and how each year was grouped in (D)–(F): before (black), between (blue), or during (light blue) marine heatwaves (MHWs).

Close modal

Phenology

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 χ2 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 and growth

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 χ2 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 (La) is:

La=Lc+(OaOc)×(LcL0)(OcO0)
1

where Lc is the length at capture, Oa is the otolith radius at age a, Oc is the otolith radius at capture, and L0 and O0 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−1) as a function of MHW designation, age, and temperature:

ln(G)HW+te(A,T,by=HW)+(A|ID)+(1|Year)
2

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−1 d−1), 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).

Phenology

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.

Table 2.

Statistical results evaluating the effect of marine heatwaves on hatch date phenology and size-at-age for Pacific cod

Response VariableFixed EffectEstimateaStd Errorat-Valueaχ2, dfbp-Valueb
Hatch date Interceptc 99.2 4.1 24.0 d – 
Betweenc −25.9 8.3 3.1 – – 
Duringc −13.9 7.2 1.9 – – 
MHWc – – – 10.9, 2 0.004 
SLe Interceptc 46.8 2.0 23.5 – – 
Betweenc −0.5 4.0 –0.1 – – 
Duringc 4.2 3.5 1.2 – – 
Agef 7.0 0.3 23.5 – – 
Between × agef 7.1 0.6 11.6 – – 
During × agef 1.7 0.7 2.4 – – 
MHWc – – – 1.7, 2 0.421 
Agef – – – 551.9, 1 <0.001 
MHW × agef – – – 134.1, 2 <0.001 
Response VariableFixed EffectEstimateaStd Errorat-Valueaχ2, dfbp-Valueb
Hatch date Interceptc 99.2 4.1 24.0 d – 
Betweenc −25.9 8.3 3.1 – – 
Duringc −13.9 7.2 1.9 – – 
MHWc – – – 10.9, 2 0.004 
SLe Interceptc 46.8 2.0 23.5 – – 
Betweenc −0.5 4.0 –0.1 – – 
Duringc 4.2 3.5 1.2 – – 
Agef 7.0 0.3 23.5 – – 
Between × agef 7.1 0.6 11.6 – – 
During × agef 1.7 0.7 2.4 – – 
MHWc – – – 1.7, 2 0.421 
Agef – – – 551.9, 1 <0.001 
MHW × agef – – – 134.1, 2 <0.001 

aResults from a linear mixed model including fixed effect estimates, standard errors, and t-values.

bResults from type-II Wald χ2 (hatch date) or type-III Wald χ2 (size-at-age) tests including χ2, degrees of freedom, and p-values.

cMarine heatwave (MHW) category was a fixed effect in both models, with the intercept representing the estimated mean value for the before MHW level, the between and during fixed effects indicating the difference between those levels and the before level, and MHW representing the fixed effect in the Wald χ2 test.

dTable cells with no value (–) distinguish the information provided between the linear mixed model and Wald χ2 tests.

eStandard length-at-capture (mm) or size-at-age.

fAge (scaled and centered) and its interaction with MHW category was a fixed effect for the size-at-age model.

Table 3.

Predicted mean values for hatch dates and standard length of Pacific cod from linear mixed models

MHW CategoryPredicted HDaPredicted SL (mm)bAge (days)
Before DOY 99 (Apr 9) c – 
During DOY 85 (Mar 26) – – 
Between DOY 73 (Mar 14) – – 
Before – 46.8d – 
During – 51.1d – 
Between – 46.4d – 
Before – 39.8e 89 
During – 42.3e 89 
Between – 32.2e 89 
Before – 44.6e 100 
During – 48.3e 100 
Between – 41.9e 100 
Before – 58.5e 132 
During – 65.6e 132 
Between – 69.9e 132 
MHW CategoryPredicted HDaPredicted SL (mm)bAge (days)
Before DOY 99 (Apr 9) c – 
During DOY 85 (Mar 26) – – 
Between DOY 73 (Mar 14) – – 
Before – 46.8d – 
During – 51.1d – 
Between – 46.4d – 
Before – 39.8e 89 
During – 42.3e 89 
Between – 32.2e 89 
Before – 44.6e 100 
During – 48.3e 100 
Between – 41.9e 100 
Before – 58.5e 132 
During – 65.6e 132 
Between – 69.9e 132 

aPredicted mean values for hatch dates (HD) as day of year (DOY) for individuals captured in years before, during, or between marine heatwaves (MHWs) from a linear mixed model that incorporated year as a random intercept.

bPredicted mean values for standard length (SL) for individuals captured in years before, during, or between MHWs from a linear mixed model that incorporated year as a random intercept and the interaction between MHW and age.

cTable cells with no value (–) distinguish the information incorporated by each model.

dPredicted mean SL values when the effect of age is ignored (main effect of MWH).

ePredicted mean SL values for three ages before, during, and between MHWs representing the youngest (89 days) and oldest (132 days) ages that were captured in both before and since periods, as well as the oldest age examined in the growth analysis (100 days).

Figure 3.

Hatch phenology of age-0 Pacific cod in the Gulf of Alaska before, during, and between marine heatwaves. (A) Individual hatch dates presented as box-and-whisker plots, where points are raw data, the box is the interquartile range, the horizontal line is the median, the vertical “whiskers” indicate the 95% confidence interval, and box colors indicate if years were before (dark gray), during (light blue), or between (blue) marine heatwaves (MHWs). (B) Density plot indicating the hatch dates before (dark gray), during (light blue), and between (blue) MHWs.

Figure 3.

Hatch phenology of age-0 Pacific cod in the Gulf of Alaska before, during, and between marine heatwaves. (A) Individual hatch dates presented as box-and-whisker plots, where points are raw data, the box is the interquartile range, the horizontal line is the median, the vertical “whiskers” indicate the 95% confidence interval, and box colors indicate if years were before (dark gray), during (light blue), or between (blue) marine heatwaves (MHWs). (B) Density plot indicating the hatch dates before (dark gray), during (light blue), and between (blue) MHWs.

Close modal

Size-at-age

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.

Figure 4.

Standard length of age-0 Pacific cod in the Gulf of Alaska before, during, and between marine heatwaves. (A) Individual standard length (SL, mm) at capture presented as box-and-whisker plots, where points are raw data, the box is the interquartile range, the horizontal line is the median, the vertical “whiskers” indicate the 95% confidence interval, and box colors indicate if years were before (dark gray), during (light blue), or between (blue) marine heatwaves (MHWs). (B) The predicted relationship from a linear mixed model between age and SL at capture within the periods before (black), during (light blue), and between (blue) MHWs. Points are raw data, and shaded regions are 95% confidence intervals.

Figure 4.

Standard length of age-0 Pacific cod in the Gulf of Alaska before, during, and between marine heatwaves. (A) Individual standard length (SL, mm) at capture presented as box-and-whisker plots, where points are raw data, the box is the interquartile range, the horizontal line is the median, the vertical “whiskers” indicate the 95% confidence interval, and box colors indicate if years were before (dark gray), during (light blue), or between (blue) marine heatwaves (MHWs). (B) The predicted relationship from a linear mixed model between age and SL at capture within the periods before (black), during (light blue), and between (blue) MHWs. Points are raw data, and shaded regions are 95% confidence intervals.

Close modal

Growth-at-age

Growth-at-age responded to temperature differently among the MHW periods (GAMM results: adjusted R2 = 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).

Table 4.

GAMM results for the response of Pacific cod growth to temperature and age before, during, and between marine heatwaves

Predictor VariableedfaFp-Value
f(age, temperature): beforeb 19.47 414.0 <0.001 
f(age, temperature): during 16.18 348.5 <0.001 
f(age, temperature): between 19.09 142.8 <0.001 
Predictor VariableedfaFp-Value
f(age, temperature): beforeb 19.47 414.0 <0.001 
f(age, temperature): during 16.18 348.5 <0.001 
f(age, temperature): between 19.09 142.8 <0.001 

aEffective degrees of freedom.

bInteractions between temperature, age, and marine heatwave, included as tensor product smooths. All smooths used cubic regression splines. GAMMs also included the random intercept of fish ID with a random slope of Age and the random intercept of Year (see Methods).

Figure 5.

GAMM-predicted relationships between age-0 Pacific cod back-calculated growth and temperature influenced by age. The general additive mixed model (GAMM) relationship presented for (A) before, (B) during, and (C) between marine heatwaves (MHWs). In the left panels, color contours indicate the Pacific cod growth rate (mm d−1) at temperature and age, with lighter colors indicating higher growth rates. In the right panels, semi-transparent black points are the data used in the model. White space outside of the contours indicates that those temperatures and ages were far outside the data (±1.25°C experienced at age).

Figure 5.

GAMM-predicted relationships between age-0 Pacific cod back-calculated growth and temperature influenced by age. The general additive mixed model (GAMM) relationship presented for (A) before, (B) during, and (C) between marine heatwaves (MHWs). In the left panels, color contours indicate the Pacific cod growth rate (mm d−1) at temperature and age, with lighter colors indicating higher growth rates. In the right panels, semi-transparent black points are the data used in the model. White space outside of the contours indicates that those temperatures and ages were far outside the data (±1.25°C experienced at age).

Close modal
Figure 6.

Absolute and relative growth of Pacific cod juveniles. (A) Daily absolute growth (mm d−1) at age and (B) relative growth (mm mm−1 d−1) at age, with lines representing annual loess smooths of collected data. (C–H) GAMM-predicted relationships between temperature and absolute growth (C, E, G) or relative growth (E, F, H) at day 10 (C, D), 40 (E, F), and 75 (G, H). In all panels, line colors represent before (black), during (light blue), and between (blue) marine heatwaves (MHWs). Gray bands around the lines represent 95% confidence intervals. For panels (C–H), temperatures where confidence intervals do not overlap indicate significant differences in growth.

Figure 6.

Absolute and relative growth of Pacific cod juveniles. (A) Daily absolute growth (mm d−1) at age and (B) relative growth (mm mm−1 d−1) at age, with lines representing annual loess smooths of collected data. (C–H) GAMM-predicted relationships between temperature and absolute growth (C, E, G) or relative growth (E, F, H) at day 10 (C, D), 40 (E, F), and 75 (G, H). In all panels, line colors represent before (black), during (light blue), and between (blue) marine heatwaves (MHWs). Gray bands around the lines represent 95% confidence intervals. For panels (C–H), temperatures where confidence intervals do not overlap indicate significant differences in growth.

Close modal

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−1 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−1 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−1and 0.46 mm d−1 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.

Hatching phenology

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 and growth

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 δ18O 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 and larger implications

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.

Abookire
,
AA
,
Litzow
,
MA
,
Malick
,
MJ
,
Laurel
,
BJ.
2022
.
Post-settlement abundance, condition, and survival in a climate-stressed population of Pacific cod
.
Canadian Journal of Fisheries and Aquatic Sciences
79
:
958
968
. DOI: http://dx.doi.org/10.1139/cjfas-2021-0224.
Almeida
,
LZ
,
Grayson
,
J
,
Ludsin
,
SA
,
Dabrowski
,
K
,
Marschall
,
EA.
2023
.
Experiential legacies of early-life dietary polyunsaturated fatty acid content on juvenile Walleye: Potential impacts from climate change
.
Ecology of Freshwater Fish
32
(
1
):
23
36
. DOI: http://dx.doi.org/10.1111/eff.12666.
Almeida
,
LZ
,
Hovick
,
SM
,
Ludsin
,
SA
,
Marschall
,
EA.
2021
.
Which factors determine the long-term effect of poor early-life nutrition? A meta-analytic review
.
Ecosphere
12
(
8
):
e03694
. DOI: http://dx.doi.org/10.1002/ecs2.3694.
Arimitsu
,
ML
,
Piatt
,
JF
,
Hatch
,
S
,
Suryan
,
RM
,
Batten
,
S
,
Bishop
,
MA
,
Campbell
,
RW
,
Coletti
,
H
,
Cushing
,
D
,
Gorman
,
K
,
Hopcroft
,
RR
,
Kuletz
,
KJ
,
Marsteller
,
C
,
McKinstry
,
C
,
McGowan
,
D
,
Moran
,
J
,
Pegau
,
S
,
Schaefer
,
A
,
Schoen
,
S
,
Straley
,
J
,
von Biela
,
VR.
2021
.
Heatwave-induced synchrony within forage fish portfolio disrupts energy flow to top pelagic predators
.
Global Change Biology
27
(
9
):
1859
1878
. DOI: http://dx.doi.org/10.1111/gcb.15556.
Atkinson
,
D.
1994
.
Temperature and organism size: A biological law for ectotherms?
Advances in Ecological Research
25
:
1
58
.
Audzijonyte
,
A
,
Fulton
,
E
,
Haddon
,
M
,
Helidoniotis
,
F
,
Hobday
,
AJ
,
Kuparinen
,
A
,
Morrongiello
,
J
,
Smith
,
ADM
,
Upston
,
J
,
Waples
,
RS.
2016
.
Trends and management implications of human-influenced life-history changes in marine ectotherms
.
Fish and Fisheries
17
(
4
):
1005
1028
. DOI: http://dx.doi.org/10.1111/faf.12156.
Audzijonyte
,
A
,
Jakubavičiūtė
,
E
,
Lindmark
,
M
,
Richards
,
SA.
2022
.
Mechanistic temperature-size rule explanation should reconcile physiological and mortality responses to temperature
.
The Biological Bulletin
243
(
2
):
220
238
. DOI: http://dx.doi.org/10.1086/722027.
Barbeaux
,
SJ
,
Ferriss
,
B
,
Palsson
,
W
,
Shotwell
,
K
,
Spies
,
I
,
Wang
,
M
,
Zandor
,
S.
2020a
.
Assessment of the Pacific cod stock in the Gulf of Alaska. North Pacific Fishery Management Council Gulf of Alaska stock assessments and fishery evaluation reports
.
Available at
https://www.fisheries.noaa.gov/resource/data/2020-assessment-pacific-cod-stock-gulf-alaska.
Accessed November 11, 2021
.
Barbeaux
,
SJ
,
Holsman
,
K
,
Zador
,
S.
2020
b.
Marine heatwave stress test of ecosystem-based fisheries management in the Gulf of Alaska Pacific cod fishery
.
Frontiers in Marine Science
7
:
703
. DOI: http://dx.doi.org/10.3389/fmars.2020.00703.
Bates
,
D
,
Mächler
,
M
,
Bolker
,
BM
,
Walker
,
SC.
2015
.
Fitting linear mixed-effects models using lme4
.
Journal of Statistical Software
67
(
1
):
1
48
. DOI: http://dx.doi.org/10.18637/jss.v067.i01.
Berkeley
,
SA
,
Chapman
,
C
,
Sogard
,
SM.
2004
.
Maternal age as a determinant of larval growth and survival in a marine fish, Sebastes melanops
.
Ecology
85
(
5
):
1258
1264
. DOI: http://dx.doi.org/10.1890/03-0706.
Bian
,
X
,
Zhang
,
X
,
Sakurai
,
Y
,
Jin
,
X
,
Gao
,
T
,
Wan
,
R
,
Yamamoto
,
J.
2014
.
Envelope surface ultrastructure and specific gravity of artificially fertilized Pacific cod Gadus macrocephalus eggs
.
Journal of Fish Biology
84
(
2
):
403
421
. DOI: http://dx.doi.org/10.1111/jfb.12292.
Bond
,
NA
,
Cronin
,
MF
,
Freeland
,
H
,
Mantua
,
N.
2015
.
Causes and impacts of the 2014 warm anomaly in the NE Pacific
.
Geophysical Research Letters
42
(
9
):
3414
3420
. DOI: http://dx.doi.org/10.1002/2015GL063306.
Bruno
,
JF
,
Carr
,
LA
,
O’Connor
,
MI.
2015
.
Exploring the role of temperature in the ocean through metabolic scaling
.
Ecology
96
(
12
):
3126
3140
. DOI: http://dx.doi.org/10.1890/14-1954.1.
Burt
,
JM
,
Hinch
,
SG
,
Patterson
,
DA.
2011
.
The importance of parentage in assessing temperature effects on fish early life history: A review of the experimental literature
.
Reviews in Fish Biology and Fisheries
21
:
377
406
. DOI: http://dx.doi.org/10.1007/s11160-010-9179-1.
Burton
,
T
,
Metcalfe
,
NB
.
2014
.
Can environmental conditions experienced in early life influence future generations?
Proceedings of the Royal Society B: Biological Sciences
281
(
1785
):
20140311
.
Campana
,
SE.
1990
.
How reliable are growth back-calculations based on otoliths?
Canadian Journal of Fisheries and Aquatic Science
47
:
2219
2227
.
Cannaday
,
JD
,
Farmer
,
TM.
2022
.
Assessing the effects of sub-tropical winter thermal conditions on coolwater fish reproduction
.
Ecology of Freshwater Fish
31
(
2
):
300
316
. DOI: http://dx.doi.org/10.1111/eff.12631.
Caputi
,
N
,
Kangas
,
MI
,
Chandrapavan
,
A
,
Hart
,
A
,
Feng
,
M
,
Marin
,
M
,
de Lestang
,
S
.
2019
.
Factors affecting the recovery of invertebrate stocks from the 2011 Western Australian extreme marine heatwave
.
Frontiers in Marine Science
6
:
484
. DOI: http://dx.doi.org/10.3389/fmars.2019.00484.
Chelgren
,
ND
,
Rosenberg
,
DK
,
Heppell
,
SS
,
Gitelman
,
AI.
2006
.
Carryover aquatic effects on survival of metamorphic frogs during pond emigration
.
Ecological Applications
16
(
1
):
250
261
. DOI: http://dx.doi.org/10.1890/04-0329.
Ciannelli
,
L
,
Neuheimer
,
AB
,
Stige
,
LC
,
Frank
,
KT
,
Durant
,
JM
,
Hunsicker
,
M
,
Rogers
,
LA
,
Porter
,
S
,
Ottersen
,
G
,
Yaragina
,
NA.
2022
.
Ontogenetic spatial constraints of sub-arctic marine fish species
.
Fish and Fisheries
23
(
2
):
342
357
. DOI: http://dx.doi.org/10.1111/faf.12619.
Croisetière
,
S
,
Bernatchez
,
L
,
Belhumeur
,
P.
2010
.
Temperature and length-dependent modulation of the MH class IIβ gene expression in brook charr (Salvelinus fontinalis) by a cis-acting minisatellite
.
Molecular Immunology
47
(
9
):
1817
1829
. DOI: http://dx.doi.org/10.1016/j.molimm.2009.12.012.
Crozier
,
LG
,
Hutchings
,
JA.
2014
.
Plastic and evolutionary responses to climate change in fish
.
Evolutionary Applications
7
(
1
):
68
87
. DOI: https://dx.doi.org/10.1111/eva.12135.
Cunningham
,
KM
,
Canino
,
MF
,
Spies
,
IB
,
Hauser
,
L.
2009
.
Genetic isolation by distance and localized fjord population structure in Pacific cod (Gadus macrocephalus): Limited effective dispersal in the northeastern Pacific Ocean
.
Canadian Journal of Fisheries and Aquatic Sciences
66
(
1
):
153
166
. DOI: http://dx.doi.org/10.1139/F08-199.
Cushing
,
DH.
1969
.
The regularity of the spawning season of some fishes
.
ICES Journal of Marine Science
33
(
1
):
81
92
. DOI: http://dx.doi.org/10.1093/icesjms/33.1.81.
Cushing
,
DH.
1990
.
Plankton production and year-class strength in fish populations: An update of the match/mismatch hypothesis
.
Advances in Marine Biology
26
:
249
293
.
Dahlke
,
FT
,
Leo
,
E
,
Mark
,
FC
,
Pörtner
,
HO
,
Bickmeyer
,
U
,
Frickenhaus
,
S
,
Storch
,
D.
2017
.
Effects of ocean acidification increase embryonic sensitivity to thermal extremes in Atlantic cod, Gadus morhua
.
Global Change Biology
23
(
4
):
1499
1510
. DOI: http://dx.doi.org/10.1111/gcb.13527.
Deutsch
,
C
,
Penn
,
JL
,
Verberk
,
WCEP
,
Inomura
,
K
,
Endress
,
M-G
,
Payne
,
JL.
2022
.
Impact of warming on aquatic body sizes explained by metabolic scaling from microbes to macrofauna
.
Proceedings of the National Academy of the Sciences
119
(
28
):
e2201345119
. DOI: http://doi.org/10.1073/pnas.2201345119.
Dunn
,
JR
,
Matarese
,
AC.
1987
.
A review of the early life history of northeast Pacific gadoid fishes
.
Fisheries Research
5
(
2–3
):
163
184
.
Durant
,
JM
,
Hjermann
,
D
,
Ottersen
,
G
,
Stenseth
,
NC.
2007
.
Climate and the match or mismatch between predator requirements and resource availability
.
Climate Research
33
(
3
):
271
283
. DOI: http://dx.doi.org/10.3354/cr033271.
Edwards
,
M
,
Richardson
,
AJ.
2004
.
Impact of climate change on marine pelagic phenology and trophic mismatch
.
Nature
430
:
881
884
. DOI: http://dx.doi.org/10.1038/nature02808.
Fedewa
,
EJ
,
Miller
,
JA
,
Hurst
,
TP
,
Jiang
,
D.
2017
.
The potential effects of pre-settlement processes on post-settlement growth and survival of juvenile northern rock sole (Lepidopsetta polyxystra) in Gulf of Alaska nursery habitats
.
Estuarine, Coastal and Shelf Science
189
:
46
57
. DOI: http://dx.doi.org/10.1016/j.ecss.2017.02.028.
Fiksen
,
Ø
,
Reglero
,
P.
2022
.
Atlantic bluefin tuna spawn early to avoid metabolic meltdown in larvae
.
Ecology
103
(
1
). DOI: http://dx.doi.org/10.1002/ecy.3568.
Fox
,
J
,
Weisberg
,
S.
2019
.
An {R} Companion to Applied Regression.
Available at
https://socialsciences.mcmaster.ca/jfox/Books/Companion/.
Accessed February 24, 2022
.
Fuiman
,
LA
,
Perez
,
KO.
2015
.
Metabolic programming mediated by an essential fatty acid alters body composition and survival skills of a marine fish
.
Proceedings of the Royal Society B: Biological Sciences
282
(
1819
):
58
64
. DOI: http://dx.doi.org/10.1098/rspb.2015.1414.
Hector
,
KL
,
Nakagawa
,
S.
2012
.
Quantitative analysis of compensatory and catch-up growth in diverse taxa
.
Journal of Animal Ecology
81
(
3
):
583
593
. DOI: http://dx.doi.org/10.1111/j.1365-2656.2011.01942.x.
Herfindal
,
I
,
van de Pol
,
M
,
Nielsen
,
JT
,
Sæther
,
B-E
,
Møller
,
AP.
2015
.
Climatic conditions cause complex patterns of covariation between demographic traits in a long-lived raptor
.
Journal of Animal Ecology
84
(
3
):
702
711
. DOI: http://dx.doi.org/10.1111/1365-2656.12318.
Hinckley
,
S
,
Stockhausen
,
WT
,
Coyle
,
KO
,
Laurel
,
BJ
,
Gibson
,
GA
,
Parada
,
C
,
Hermann
,
AJ
,
Doyle
,
MJ
,
Hurst
,
TP
,
Punt
,
AE
,
Ladd
,
C.
2019
.
Connectivity between spawning and nursery areas for Pacific cod (Gadus macrocephalus) in the Gulf of Alaska
.
Deep-Sea Research Part II: Topical Studies in Oceanography
165
(
April
):
113
126
. DOI: http://dx.doi.org/10.1016/j.dsr2.2019.05.007.
Hobday
,
AJ
,
Alexander
,
LV
,
Perkins
,
SE
,
Smale
,
DA
,
Straub
,
SC
,
Oliver
,
ECJ
,
Benthuysen
,
JA
,
Burrows
,
MT
,
Donat
,
MG
,
Feng
,
M
,
Holbrook
,
NJ
,
Moore
,
PJ
,
Scannell
,
HA
,
Gupta
,
A
,
Wernberg
,
T.
2016
.
A hierarchical approach to defining marine heatwaves
.
Progress in Oceanography
141
:
227
238
. DOI: http://dx.doi.org/10.1016/j.pocean.2015.12.014.
Hobday
,
AJ
,
Oliver
,
ECJ
,
Sen Gupta
,
A
,
Benthuysen
,
JA
,
Burrows
,
MT
,
Donat
,
MG
,
Holbrook
,
NJ
,
Moore
,
PJ
,
Thomsen
,
MS
,
Wernberg
,
T
,
Smale
,
DA.
2018
.
Categorizing and naming marine heatwaves
.
Oceanography
31
(
2
):
162
173
.
Hurst
,
TP
,
Cooper
,
DW
,
Scheingross
,
JS
,
Seale
,
EM
,
Laurel
,
BJ
,
Spencer
,
ML.
2009
.
Effects of ontogeny, temperature, and light on vertical movements of larval Pacific cod (Gadus macrocephalus)
.
Fisheries Oceanography
18
(
5
):
301
311
. DOI: http://dx.doi.org/10.1111/j.1365-2419.2009.00512.x.
Hurst
,
TP
,
Laurel
,
BJ
,
Ciannelli
,
L.
2010
.
Ontogenetic patterns and temperature-dependent growth rates in early life stages of Pacific cod (Gadus macrocephalus)
.
Fishery Bulletin
108
(
4
):
382
392
.
Hurst
,
TP
,
Munch
,
SB
,
Lavelle
,
KA.
2012
.
Thermal reaction norms for growth vary among cohorts of Pacific cod (Gadus macrocephalus)
.
Marine Biology
159
(
10
):
2173
2183
. DOI: http://dx.doi.org/10.1007/s00227-012-2003-9.
Jobling
,
M.
1997
.
Temperature and growth: Modulation of growth rate via temperature change
, in
Wood
,
CM
,
McDonald
,
DG
eds.,
Global warming: Implications for freshwater and marine fish
.
Cambridge, UK
:
Cambridge University Press
:
225
254
. DOI: http://dx.doi.org/10.1017/cbo9780511983375.010.
Jørgensen
,
KEM
,
Neuheimer
,
AB
,
Jorde
,
PE
,
Knutsen
,
H
,
Grønkjær
,
P.
2020
.
Settlement processes induce differences in daily growth rates between two co-existing ecotypes of juvenile cod Gadus morhua
.
Marine Ecology Progress Series
650
:
175
189
. DOI: http://dx.doi.org/10.3354/meps13433.
Jutfelt
,
F
,
Norin
,
T
,
Åsheim
,
ER
,
Rowsey
,
LE
,
Andreassen
,
AH
,
Morgan
,
R
,
Clark
,
TD
,
Speers-Roesch
,
B.
2021
.
‘Aerobic scope protection’ reduces ectotherm growth under warming
.
Functional Ecology
35
(
7
):
1397
1407
. DOI: http://dx.doi.org/10.1111/1365-2435.13811.
Kastelle
,
CR
,
Helser
,
TE
,
Laurel
,
BJ
,
Copeman
,
LA
,
Stone
,
KR
,
McKay
,
JL.
2022
.
Oxygen isotope fractionation in otoliths: Experimental results from four North Pacific and Arctic gadid species
.
Marine Ecology Progress Series
686
:
159
175
. DOI: http://dx.doi.org/10.3354/meps13985.
Ketchen
,
KS.
1961
.
Observations on the ecology of the Pacific cod (Gadus macrocephalus) in Canadian waters
.
Journal of the Fisheries Research Board of Canada
18
(
4
):
513
558
.
Kjesbu
,
OS
,
Righton
,
D
,
Krüger-Johnsen
,
M
,
Thorsen
,
A
,
Michalsen
,
K
,
Fonn
,
M
,
Witthames
,
PR.
2010
.
Thermal dynamics of ovarian maturation in Atlantic cod (Gadus morhua)
.
Canadian Journal of Fisheries and Aquatic Sciences
67
(
4
):
605
625
. DOI: http://dx.doi.org/10.1139/f10-011.
Koussoroplis
,
AM
,
Nussbaumer
,
J
,
Arts
,
MT
,
Guschina
,
IA
,
Kainz
,
MJ.
2014
.
Famine and feast in a common freshwater calanoid: Effects of diet and temperature on fatty acid dynamics of Eudiaptomus gracilis
.
Limnology and Oceanography
59
(
3
):
947
958
. DOI: http://dx.doi.org/10.4319/lo.2014.59.3.0947.
Laufkötter
,
C
,
Zscheischler
,
J
,
Frölicher
,
TL.
2020
.
High-impact marine heatwaves attributable to human-induced global warming
.
Science
369
(
6511
):
1621
1625
.
Laurel
,
BJ
,
Abookire
,
AA
,
Barbeaux
,
SJ
,
Almeida
,
LZ
,
Copeman
,
LA
,
Duffy-Anderson
,
J
,
Hurst
,
TP
,
Litzow
,
MA
,
Kristiansen
,
T
,
Miller
,
JA
,
Palsson
,
W
,
Rooney
,
S
,
Thalmann
,
HL
,
Rogers
,
LA
.
2023
.
Pacific cod in the Anthropocene: An early life history perspective under changing thermal habitats
.
Fish and Fisheries
24
(
6
):
959
978
. DOI: http://dx.doi.org/10.1111/faf.12779.
Laurel
,
BJ
,
Cote
,
D
,
Gregory
,
RS
,
Rogers
,
L
,
Knutsen
,
H
,
Olsen
,
EM.
2017
.
Recruitment signals in juvenile cod surveys depend on thermal growth conditions
.
Canadian Journal of Fisheries and Aquatic Sciences
74
(
4
):
511
523
. DOI: http://dx.doi.org/10.1139/cjfas-2016-0035.
Laurel
,
BJ
,
Hunsicker
,
ME
,
Ciannelli
,
L
,
Hurst
,
TP
,
Duffy-Anderson
,
J
,
O’Malley
,
R
,
Behrenfeld
,
M.
2021
.
Regional warming exacerbates match/mismatch vulnerability for cod larvae in Alaska
.
Progress in Oceanography
193
: 102555. DOI: http://dx.doi.org/10.1016/j.pocean.2021.102555.
Laurel
,
BJ
,
Hurst
,
TP
,
Ciannelli
,
L.
2011
.
An experimental examination of temperature interactions in the match–mismatch hypothesis for Pacific cod larvae
.
Canadian Journal of Fisheries and Aquatic Sciences
68
(
1
):
51
61
. DOI: http://dx.doi.org/10.1139/F10-130.
Laurel
,
BJ
,
Hurst
,
TP
,
Copeman
,
LA
,
Davis
,
MW.
2008
.
The role of temperature on the growth and survival of early and late hatching Pacific cod larvae (Gadus macrocephalus)
.
Journal of Plankton Research
30
(
9
):
1051
1060
. DOI: http://dx.doi.org/10.1093/plankt/fbn057.
Laurel
,
BJ
,
Knoth
,
BA
,
Ryer
,
CH.
2016
.
Growth, mortality, and recruitment signals in age-0 gadids settling in coastal Gulf of Alaska
.
ICES Journal of Marine Science
73
(
9
):
2227
2237
.
Laurel
,
BJ
,
Rogers
,
LA.
2020
.
Loss of spawning habitat and prerecruits of Pacific cod during a Gulf of Alaska heatwave
.
Canadian Journal of Fisheries and Aquatic Sciences
77
(
4
):
644
650
. DOI: http://dx.doi.org/10.1139/cjfas-2019-0238.
Laurel
,
BJ
,
Stoner
,
AW
,
Ryer
,
CH
,
Hurst
,
TP
,
Abookire
,
AA.
2007
.
Comparative habitat associations in juvenile Pacific cod and other gadids using seines, baited cameras and laboratory techniques
.
Journal of Experimental Marine Biology and Ecology
351
(
1–2
):
42
55
. DOI: http://dx.doi.org/10.1016/j.jembe.2007.06.005.
Lenth
,
RV.
2021
.
Emmeans: Estimated marginal means, aka least-squares means
.
Available at
https://cran.r-project.org/package=emmeans.
Accessed February 24, 2022
.
Li
,
Z
,
Yamamoto
,
J
,
Sakurai
,
Y.
2015
.
Vertical position, specific gravity and swimming ability of Pacific cod Gadus macrocephalus yolk-sac larvae reared at four temperatures
.
Fisheries Science
81
(
5
):
883
889
. DOI: http://dx.doi.org/10.1007/s12562-015-0911-6.
Litz
,
MNC
,
Miller
,
JA
,
Copeman
,
LA
,
Hurst
,
TP.
2017
.
Effects of dietary fatty acids on juvenile salmon growth, biochemistry, and aerobic performance: A laboratory rearing experiment
.
Journal of Experimental Marine Biology and Ecology
494
:
20
31
. DOI: http://dx.doi.org/10.1016/j.jembe.2017.04.007.
Litzow
,
MA
,
Abookire
,
AA
,
Duffy-Anderson
,
JT
,
Laurel
,
BJ
,
Malick
,
MJ
,
Rogers
,
LA.
2022
.
Predicting year class strength for climate-stressed gadid stocks in the Gulf of Alaska
.
Fisheries Research
249
(
S39–S43
):
106250
. DOI: http://dx.doi.org/10.1016/j.fishres.2022.106250.
Litzow
,
MA
,
Hunsicker
,
ME
,
Ward
,
EJ
,
Anderson
,
SC
,
Gao
,
J
,
Zador
,
SG
,
Batten
,
S
,
Dressel
,
SC
,
Duffy-Anderson
,
J
,
Fergusson
,
E
,
Hopcroft
,
RR
,
Laurel
,
BJ
,
O’Malley
,
R.
2020
.
Evaluating ecosystem change as Gulf of Alaska temperature exceeds the limits of preindustrial variability
.
Progress in Oceanography
186
(
February
):
102393
. DOI: http://dx.doi.org/10.1016/j.pocean.2020.102393.
Litzow
,
MA
,
Malick
,
MJ
,
Abookire
,
AA
,
Duffy-Anderson
,
J
,
Laurel
,
BJ
,
Ressler
,
PH
,
Rogers
,
LA.
2021
.
Using a climate attribution statistic to inform judgments about changing fisheries sustainability
.
Scientific Reports
11
(
1
):
23924
. DOI: http://dx.doi.org/10.1038/s41598-021-03405-6.
Meffe
,
GK
,
Weeks
,
SC
,
Mulvey
,
M
,
Kandl
,
KL.
1995
.
Genetic differences in thermal tolerance of eastern mosquitofish (Gambusia holbrooki; Poeciliidae) from ambient and thermal ponds
.
Canadian Journal of Fisheries and Aquatic Sciences
52
:
2704
2711
.
Menzel
,
A
,
Sparks
,
TH
,
Estrella
,
N
,
Koch
,
E
,
Aasas
,
A
,
Ahass
,
R
,
Alm-Kubler
,
K
,
Bissolli
,
P
,
Braslavska
,
O
,
Briede
,
A
,
Chmielewski
,
FM
,
Črepinšek
,
Z
,
Curnel
,
Y
,
Dahl
,
A
,
Defila
,
C
,
Donnelly
,
A
,
Filella
,
I
,
Jatczak
,
K
,
Mage
,
F
,
Mestre
,
A
,
Nordli
,
O
,
Penuelas
,
J
,
Pirinen
,
P
,
Remisova
,
V
,
Scheifinger
,
H
,
Striz
,
M
,
Susnik
,
A
,
Emil Wielgolaski
,
F
,
Zach
,
S
,
Žust
,
A
.
2006
.
European phenological response to climate change matches the warming pattern
.
Global Change Biology
12
:
1969
1976
. DOI: http://dx.doi.org/10.1111/j.1365-2486.2006.01193.x.
Miller
,
TJ
,
Crowder
,
LB
,
Rice
,
JA
,
Marschall
,
EA.
1988
.
Larval size and recruitment mechanisms in fishes: Toward a conceptual framework
.
Canadian Journal of Fisheries and Aquatic Sciences
45
(
9
):
1657
1670
. DOI: http://dx.doi.org/10.1139/f88-197.
Morrongiello
,
JR
,
Thresher
,
RE.
2015
.
A statistical framework to explore ontogenetic growth variation among individuals and populations: A marine fish example
.
Ecological Monographs
85
(
1
):
93
115
. DOI: http://dx.doi.org/10.1890/13-2355.1.
Morrongiello
,
JR
,
Walsh
,
CT
,
Gray
,
CA
,
Stocks
,
JR
,
Crook
,
DA.
2014
.
Environmental change drives long-term recruitment and growth variation in an estuarine fish
.
Global Change Biology
20
(
6
):
1844
1860
. DOI: http://dx.doi.org/10.1111/gcb.12545.
Murphy
,
HM
,
Warren-Myers
,
FW
,
Jenkins
,
GP
,
Hamer
,
PA
,
Swearer
,
SE.
2014
.
Variability in size-selective mortality obscures the importance of larval traits to recruitment success in a temperate marine fish
.
Oecologia
175
(
4
):
1201
1210
. DOI: http://dx.doi.org/10.1007/s00442-014-2968-9.
Narimatsu
,
Y
,
Hattori
,
T
,
Ueda
,
Y
,
Matsuzaka
,
H
,
Shiogaki
,
M.
2007
.
Somatic growth and otolith microstructure of larval and juvenile Pacific cod Gadus macrocephalus
.
Fisheries Science
73
(
6
):
1257
1264
. DOI: http://dx.doi.org/10.1111/j.1444-2906.2007.01463.x.
Neidetcher
,
SK
,
Hurst
,
TP
,
Ciannelli
,
L
,
Logerwell
,
EA.
2014
.
Spawning phenology and geography of Aleutian Islands and eastern Bering Sea Pacific cod (Gadus macrocephalus)
.
Deep-Sea Research Part II: Topical Studies in Oceanography
109
:
204
214
. DOI: http://dx.doi.org/10.1016/j.dsr2.2013.12.006.
Neuheimer
,
AB
,
MacKenzie
,
BR
,
Payne
,
MR.
2018
.
Temperature-dependent adaptation allows fish to meet their food across their species’ range
.
Science Advances
4
(
7
). DOI: http://dx.doi.org/10.1126/sciadv.aar4349.
Ohlberger
,
J
,
Thackeray
,
SJ
,
Winfield
,
IJ
,
Maberly
,
SC
,
Vøllestad
,
LA.
2014
.
When phenology matters: Age–size truncation alters population response to trophic mismatch
.
Proceedings of the Royal Society B: Biological Sciences
281
(
1793
). DOI: http://dx.doi.org/10.1098/rspb.2014.0938.
Paradis
,
AR
,
Pepin
,
P
,
Brown
,
JA.
1996
.
Vulnerability of fish eggs and larvae to predation: Review of the influence of the relative size of prey and predator
.
Canadian Journal of Fisheries and Aquatic Sciences
53
(
6
):
1226
1235
.
Peck
,
MA
,
Huebert
,
KB
,
Llopiz
,
JK.
2012
.
Intrinsic and extrinsic factors driving match–mismatch dynamics during the early life history of marine fishes
.
Advances in Ecological Research
47
:
177
302
. DOI: http://dx.doi.org/10.1016/B978-0-12-398315-2.00003-X.
Peterson Williams
,
MJ
,
Robbins Gisclair
,
B
,
Cerny-Chipman
,
E
,
LeVine
,
M
,
Peterson
,
T.
2021
.
The heat is on: Gulf of Alaska Pacific cod and climate-ready fisheries
.
ICES Journal of Marine Science
2016
(
2
). DOI: http://dx.doi.org/10.1093/icesjms/fsab032.
Poloczanska
,
ES
,
Burrows
,
MT
,
Brown
,
CJ
,
Molinos
,
JG
,
Halpern
,
BS
,
Hoegh-Guldberg
,
O
,
Kappel
,
CV
,
Moore
,
PJ
,
Richardson
,
AJ
,
Schoeman
,
DS
,
Sydeman
,
WJ.
2016
.
Responses of marine organisms to climate change across oceans
.
Frontiers in Marine Science
3
:
62
. DOI: http://dx.doi.org/10.3389/fmars.2016.00062.
Post
,
JR
,
Evans
,
DO.
1989
.
Size-dependent ovewinter mortality of young-of-the-year yellow perch (Perca flavescens): Laboratory, in situ enclosure, and field experiments
.
Canadian Journal of Fisheries and Aquatic Sciences
46
(
11
):
1958
1968
. DOI: http://dx.doi.org/10.1139/f89-246.
Post
,
JR
,
Parkinson
,
EA.
2001
.
Energy allocation strategy in young fish: Allometry and survival
.
Ecology
82
(
4
):
1040
1051
.
R Core Team
.
2021
.
R: A Language and environment for statistical computing
.
Available at
https://www.r-project.org/.
Accessed September 12, 2021
.
Rasmussen
,
NL
,
Van Allen
,
BG
,
Rudolf
,
VHW.
2014
.
Linking phenological shifts to species interactions through size-mediated priority effects
.
Journal of Animal Ecology
83
(
5
):
1206
1215
. DOI: http://dx.doi.org/10.1111/1365-2656.12203.
Reusch
,
TBH.
2014
.
Climate change in the oceans: Evolutionary versus phenotypically plastic responses of marine animals and plants
.
Evolutionary Applications
7
(
1
):
104
122
. DOI: http://dx.doi.org/10.1111/eva.12109.
Rijnsdorp
,
AD
,
Peck
,
MA
,
Engelhard
,
GH
,
Möllmann
,
C
,
Pinnegar
,
JK.
2009
.
Resolving the effect of climate change on fish populations
.
ICES Journal of Marine Science
66
(
September
):
1570
1583
.
Rogers
,
LA
,
Dougherty
,
AB.
2019
.
Effects of climate and demography on reproductive phenology of a harvested marine fish population
.
Global Change Biology
25
(
2
):
708
720
. DOI: http://dx.doi.org/10.1111/gcb.14483.
Sakamoto
,
T
,
Takahashi
,
M
,
Chung
,
MT
,
Rykaczewski
,
RR
,
Komatsu
,
K
,
Shirai
,
K
,
Ishimura
,
T
,
Higuchi
,
T.
2022
.
Contrasting life-history responses to climate variability in eastern and western North Pacific sardine populations
.
Nature Communications
13
(
1
):
5298
. DOI: http://dx.doi.org/10.1038/s41467-022-33019-z.
Schlegel
,
RW
,
Smit
,
AJ.
2018
.
heatwaveR: A central algorithm for the detection of heatwaves and cold-spells
.
Journal of Open Source Software
3
(
27
):
821
. DOI: http://dx.doi.org/10.21105/joss.00821.
Schneider
,
KN
,
Newman
,
RM
,
Card
,
V
,
Weisberg
,
S
,
Pereira
,
DL.
2010
.
Timing of walleye spawning as an indicator of climate change
.
Transactions of the American Fisheries Society
139
(
4
):
1198
1210
. DOI: http://dx.doi.org/10.1577/t09-129.1.
Sih
,
A
,
Ferrari
,
MCO
,
Harris
,
DJ.
2011
.
Evolution and behavioural responses to human-induced rapid environmental change
.
Evolutionary Applications
4
(
2
):
367
387
. DOI: http://dx.doi.org/10.1111/j.1752-4571.2010.00166.x.
Stark
,
JW.
2007
.
Geographic and seasonal variations in maturation and growth of female Pacific cod (Gadus macrocephalus) in the Gulf of Alaska and Bering Sea
.
Fishery Bulletin
105
(
3
):
396
407
.
Suryan
,
RM
,
Arimitsu
,
ML
,
Coletti
,
HA
,
Hopcroft
,
RR
,
Lindeberg
,
MR
,
Barbeaux
,
SJ
,
Batten
,
SD
,
Burt
,
WJ
,
Bishop
,
MA
,
Bodkin
,
JL
,
Brenner
,
R
,
Campbell
,
R
,
Cushing
,
D
,
Danielson
,
S
,
Dorn
,
MW
,
Drummond
,
B
,
Esler
,
D
,
Gelatt
,
T
,
Hanselman
,
DH
,
Hatch
,
SA
,
Haught
,
S
,
Holderied
,
K
,
Iken
,
K
,
Irons
,
DB
,
Kettle
,
AB
,
Kimmel
,
DG
,
Konar
,
B
,
Kuletz
,
KJ
,
Laurel
,
B
,
Maniscalco
,
JM
,
Matkin
,
CO
,
Mckinstry
,
CAE
,
Monson
,
DH
,
Moran
,
JR
,
Olsen
,
D
,
Palsson
,
WA
,
Pegau
,
WS
,
Piatt
,
JF
,
Rogers
,
LA
,
Rojek
,
NA
,
Schaefer
,
A
,
Spies
,
IB
,
Straley
,
J
,
Strom
,
SL
,
Sweeney
,
KL
,
Szymkowiak
,
M
,
Weitzman
,
BP
,
Yasumiishi
,
EM
,
Zador
,
SG
.
2021
.
Ecosystem response persists after a prolonged marine heatwave
.
Scientific Reports
11
:
6235
. DOI: http://dx.doi.org/10.1038/s41598-021-83818-5.
Thalmann
,
HL
,
Daly
,
EA
,
Brodeur
,
RD.
2020
.
Two anomalously warm years in the Northern California current: Impacts on early marine steelhead diet composition, morphology, and potential survival
.
Transactions of the American Fisheries Society
149
(
4
):
369
382
. DOI: http://dx.doi.org/10.1002/tafs.10244.
Van Allen
,
BG
,
Rudolf
,
VHW.
2016
.
Carryover effects drive competitive dominance in spatially structured environments
.
Proceedings of the National Academy of Sciences of the United States of America
113
(
25
):
6939
6944
. DOI: http://dx.doi.org/10.1073/pnas.1520536113.
Vigliola
,
L
,
Harmelin-Vivien
,
M
,
Meekan
,
MG.
2000
.
Comparison of techniques of back-calculation of growth and settlement marks from the otoliths of three species of Diplodus from the Mediterranean Sea
.
Canadian Journal of Fisheries and Aquatic Sciences
57
:
1291
1299
.
Vigliola
,
L
,
Meekan
,
MG.
2002
.
Size at hatching and planktonic growth determine post-settlement survivorship of a coral reef fish
.
Oecologia
131
(
1
):
89
93
. DOI: http://dx.doi.org/10.1007/s00442-001-0866-4.
Vigliola
,
L
,
Meekan
,
MG.
2009
.
The back-calculation of fish growth from otoliths
, in
Green
,
BS
,
Mapstone
,
BD
,
Carlos
,
G
,
Begg
,
GA
eds.,
Tropical fish otoliths: Information for assessment, management and ecology.
Dordrecht, the Netherlands
:
Springer
:
174
211
. (
Reviews: Methods and technologies in fish biology and fisheries; vol. 11)
. DOI: http://dx.doi.org/10.1007/978-1-4020-5775-5_6.
Weingartner
,
T
,
Eisner
,
L
,
Eckert
,
GL
,
Danielson
,
S.
2009
.
Southeast Alaska: Oceanographic habitats and linkages
.
Journal of Biogeography
36
(
3
):
387
400
. DOI: http://dx.doi.org/10.1111/j.1365-2699.2008.01994.x.
Wickham
,
H.
2016
.
Ggplot2: Elegant graphics for data analysis
.
Available at
https://ggplot2.tidyverse.org.
Accessed December 20, 2020
.
Williams
,
JW
,
Jackson
,
ST.
2007
.
Novel climates, no-analog communities, and ecological surprises
.
Frontiers in Ecology and the Environment
5
(
9
):
475
482
. DOI: http://dx.doi.org/10.1890/070037.
Wood
,
SN
.
2017
.
Generalized additive models: An introduction with R
.
Boca Raton, FL
:
CRC Press
.
Wootton
,
HF
,
Morrongiello
,
JR
,
Schmitt
,
T
,
Audzijonyte
,
A.
2022
.
Smaller adult fish size in warmer water is not explained by elevated metabolism
.
Ecology Letters
25
(
5
):
1177
1188
. DOI: http://dx.doi.org/10.1111/ele.13989.
Zambonino-Infante
,
JL
,
Claireaux
,
G
,
Ernande
,
B
,
Jolivet
,
A
,
Quazuguel
,
P
,
Sévère
,
A
,
Huelvan
,
C
,
Mazurais
,
D.
2013
.
Hypoxia tolerance of common sole juveniles depends on dietary regime and temperature at the larval stage: Evidence for environmental conditioning
.
Proceedings of the Royal Society B: Biological Sciences
280
(
1758
):
20123022
. DOI: http://dx.doi.org/10.1098/rspb.2012.3022.

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

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC-BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. See http://creativecommons.org/licenses/by/4.0/.

Supplementary data