Five decades of change in somatic growth of Pacific hake from Puget Sound and Strait of Georgia

Study species

Pacific hake (Merluccius productus) is a pelagic species distributed from Baja California, USA to northern B.C., Canada (Berger et al., 2019). We studied size-at-age patterns of Pacific hake within the Georgia Basin DPS, which includes individuals from populations in greater Puget Sound (PS) and the Strait of Georgia (SoG) (Fig. 1) for which genetic differences have been reported (Iwamoto, Ford & Gustafson, 2004; Iwamoto et al., 2015). In PS, spawning aggregations occur in Port Susan (Pedersen, 1985) and Dabob Bay (Utter, Stormont & Hodgins, 1970). In SoG, spawning aggregations occur in Stuart Channel (McFarlane & Beamish, 1985) and Saanich Inlet (Beamish, Smith & Scarsbrook, 1978). Individuals mature by age 3–4 (McFarlane & Beamish, 1985) and spawning occurs primarily from February through April, peaking in March (Goñi, 1988). Hatching takes place 4–6 days after fertilization, at which point larvae are 2–3 mm total length (Bailey, 1982b; Hollowed, 1992).

Fish collections and otolith preparation

Using pelagic trawls, juvenile and adult Pacific hake were collected from PS and SoG spawning grounds by the Washington Department of Fish and Wildlife (WDFW) (Burger et al., 2019) and Fisheries and Oceans, Canada (DFO), respectively (Fig. 1) across several years (Table 1). Fish total length (TL) was measured, and otoliths were removed and stored dry. Left sagittal otoliths were mounted on glass microscope slides using thermoplastic cement; otoliths were then polished in a sagittal plane using a grinding wheel with slurries of 600-grit silicon carbide, 5.0 alumina oxide, and 1.0 micropolish (Buehler²) until the core was visible. Images of each otolith were taken with transmitted light using a digital camera (Mediacybernetics, EvolutionMP) at 4 × magnification and analyzed using image processing software (Image Pro Plus, version 7.0).

Age, length, and growth

To estimate age and change in size over time, we examined the annual pattern of otolith growth, which consists of an opaque and translucent zone (Chilton & Beamish, 1982). Under transmitted light, the opaque zone corresponds to a summertime period when fish do most of their growing, while the translucent zone is a period of reduced winter growth (Beamish, 1979; Chilton & Beamish, 1982) (Fig. 2A). For each fish, we estimated age by counting the number of otolith translucent zones (Bailey, 1982b; Beamish, 1979) along a transect perpendicular to its longest axis on the dorsal side. Age was estimated for PS and SoG fish by WDFW in Olympia, Washington, and by DFO at the Pacific Biological Station in Nanaimo, BC, respectively. A second age was estimated by the U.S. National Marine Fisheries Service, Northwest Fisheries Science Center in Seattle, Washington. If the two ages differed for an individual, then age was estimated a third and final time by the Northwest Fisheries Science Center. Fish were assigned an age and the otolith included for statistical analysis only if two age estimates were the same.

Because somatic and otolith growth are positively correlated and consistent between both populations (R² = 0.83, n = 532, Fig. S1) we used the otolith radius at opaque and translucent zones to estimate fish length at the end of each summer and winter, respectively. Specifically, to back-calculate summer and winter body sizes we measured the otolith radius to each opaque (O_s) and translucent (O_w) zone, as well as the distance from the otolith core to the edge (i.e., otolith radius at time of capture, O_c) (Beamish, 1979) (Fig. 2A). For each individual, fish length after each summer (L_s) and winter (L_w) was estimated using the Biological Intercept equation (Campana, 2011; Campana & Jones, 1992): ${\displaystyle L_s=L_c+\frac{\left(O_s-O_c\right)x\left(L_c-L_i\right)}{\left(O_c-O_i\right)}}$

where L_i and O_i were the biological intercepts and were defined as 3.9 mm and 14.1 µm, respectively and L_c corresponds to the length at capture. The biological intercepts corresponded to the size at first feeding and were based on work by Bailey (1982a) and Butler & Nishimoto (1997). In addition, we calculated the estimated amount an individual grew (mm) each summer and winter as the difference between its length at any summer or winter from its length estimated for the previous winter or summer, respectively.

Our statistical analyses of length and growth data had two objectives. Our first objective was to identify when size differences arose between PS and SoG fish. We used a repeated measures ANOVA to compare length and growth between fish collected from PS and SoG, regardless of hatch year. Repeated measures ANOVA was appropriate because our otolith-derived estimates of length and growth represent repeated measures that were calculated for each age-season of a fish’s life (i.e., first summer, first winter, second summer, etc.) (Table 2).

The second objective was to determine the extent to which length and growth have changed across hatch years, grouped by decade, within a population. We used a repeated measures ANOVA to investigate if length and growth differed among each year/season and among hatch years. We grouped fish according to the decade in which they hatched because low sample sizes for certain hatch years precluded year-specific analysis. If significant differences were detected we performed a Tukey’s post-hoc test to determine which decades were significantly different in terms of length and growth. We used a mixed effects model to describe the repeated measures analysis using the lme function in the nlme package of RStudio (version 1.1.463) (R Core Team, 2021). ‘Individual’ was treated as a random variable in the model and we included an autocorrelation structure for lags in the time variable (i.e., year/season) via the ACF function in the nlme package.

Otolith chemistry

We chemically analyzed the same otoliths from which we estimated age and seasonal growth since ratios of certain isotopes are known to reflect the environment in which the fish resides (Elsdon et al., 2008). Prior to their chemical analysis, otoliths mounted to microscope slides were cleaned by rinsing with 95% ethanol. We used laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) at the GeoAnalytical Lab, Washington State University. Each laser ablation of an otolith obtained trace elemental concentrations from the otolith edge, which corresponds in space and time to their collection from a spawning ground, to the otolith core, which corresponds to hatching.

Data were collected as previously described in Chittaro et al. (2013). Specifically, isotope concentrations were determined using high resolution single collector inductively coupled plasma mass spectrometry (Finnigan Element2, with helium as the carrier gas) that was coupled with a laser ablation system (New Wave UP-213, frequency of 20-Hz, 30-µm spot size). Using an automated microscope stage, the laser beam was focused on the otolith, and a transect from edge to its core was ablated at a speed of 15 µm s⁻¹. Each transect was placed along the same axis used to estimate age and size, and formed an edge-to-core scan line corresponding to the entire life of the fish. Data acquisition of the LA-ICP-MS lasted 240 s, 20 s of which were designated for instrument calibration and gas background counts prior to the start of each ablation. To correct for instrument drift, we obtained a glass standard doped with trace elements from the National Institute of Standards and Technology (NIST 610). We analyzed this standard both at the beginning and end of each sample set (i.e., 16–20 otoliths). Calcium was used as an internal standard to compensate for signal variation caused by differences in the mass of ablated material. We calculated detection limits for each isotope as the average background plus three standard deviations.

From the laser ablation of each core to edge transect, six isotopes (magnesium-25, Mg; calcium-43, Ca; manganese-55, Mn; strontium-86, Sr; barium-138, Ba; and lead-208, Pb) were analyzed by ICP-MS. Counts per second were measured for each isotope along this transect, from which we calculated ratios of each isotope to calcium (e.g., Mg:Ca, mm mol⁻¹) (Fig. 2B). Next, we calculated average isotope ratios for each period of summer and winter growth. Specifically, we used the otolith radius of each opaque and translucent zone (summer and winter growth, respectively) to designate the section of the core-to-edge transect for which to calculate the average for every isotope ratio (Fig. 2). When calculating the average isotopic ratio for the first summer of growth we excluded the yolk-sac stage because this stage has been shown in other species to be a product of maternal investment, and thus is not a site-specific signature (Brophy, Jeffries & Danilowicz, 2004; Chittaro et al., 2006). The yolk-sac stage was excluded from our isotope ratios by omitting the chemical signature of the laser transect segment within 100-µm of the otolith core. We were confident that this method excluded any maternal signal since yolk-sac absorption occurs 5–7 days after hatching (Bailey, 1982b) and 30–40 day old larvae are reported to have an otolith radius of 100–150 µm (Butler & Nishimoto, 1997).

Isotopes were included for statistical analyses if they met two criteria (Chittaro et al., 2013): concentrations of the isotope were greater than the detection limit in more than 80% of otoliths analyzed, and concentrations of the isotope in NIST samples were determined with satisfactory precision (coefficient of variation < 10%).

Generalized linear model

We used a generalized linear modeling (GLM) approach to investigate the extent to which variability in somatic growth (dependent variable) was explained by seven variables: year, harbor seal (Phoca vitulina) abundance, sea surface temperature, regional climate index, and several isotope ratios. This analysis was performed separately for each of three summers of growth, corresponding to the first, second, and third year of life. We used the glm function in the stats package of RStudio (R Core Team, 2021) and specified a gamma family distribution with a log link to account for the normally distributed, but positive, growth data.

We included year to account for interannual variation in marine environmental conditions not otherwise accounted for explicitly in the model. ‘Year’ corresponds to the year in which the summer growth occurred. For example, if a fish hatched in 1990 and we are investigating its second summer of growth, then the year associated with this growth would be 1992. Harbor seal abundance (Chasco et al., 2017) was included to account for the known predation pressure of harbor seals on Pacific hake (Gustafson et al., 2000). According to Saunders & McFarlane (1999), SoG hake represented 42% of the diets in harbor seals. Because Pacific hake are ectothermic their metabolic demands are positively correlated to water temperature. We therefore included sea surface temperature, averaged from April to August, collected from a buoy located at Race Rocks on the southeastern end of Victoria Island, British Columbia, Canada. To account for regional climate patterns we included Pacific Northwest Index (PNI) (Ebbesmeyer & Strickland 1995), which is a composite of location-specific air temperature, total precipitation, and snowpack data. Finally, because some aspects of water quality are recorded in the uptake of trace elements onto the growing otolith surface we used the chemical constituents of Pacific hake otoliths as a means to gain insight into the habitat where they resided during each summer. Specifically, we used average Mn:Ca as an indicator of hypoxic conditions (Limburg & Casini, 2018; Limburg et al., 2011; Limburg et al., 2015; Mohan & Walther, 2016), average Ba:Ca as an indicator of nutrient rich upwelling (Bath et al., 2000; Mohan et al., 2018), and average Sr:Ca (Macdonald & Crook, 2010; Martin, Thorrold & Jones, 2004; Zimmerman, 2005) and Ba:Ca (Elsdon & Gillanders, 2005; (Hamer, Jenkins & Coutin, 2006; Miller, Gray & Merz, 2010) as an indicator of salinity. Lead was excluded from our analyses because little is understood about how it relates to water quality.

We ran GLMs only for periods of summer growth and separately for each of the first three years, and each population. We limited our analyses to summer growth because yearly changes in growth were more dramatic, relative to winter, (see results) and harbor seal abundance data were collected in spring and summer (Chasco et al., 2017). The most complex model was describes as: ${\displaystyle G\sim Yr+Seal+Temp+PNI+Mn+Sr+Ba}$

where G is the summer growth estimated from otoliths, Yr is the year the summer growth occurred, Seal is harbor seal abundance, Temp is the sea surface temperature, PNI is the Pacific Northwest Index, and Mn, Sr, and Ba are the average Mn:Ca, Sr:Ca, and Ba:Ca, respectively, corresponding to the summer of otolith growth (e.g., first summer of growth).

All possible GLM model combinations (n = 127) were run and model parameters were estimated by maximizing the likelihood function. To compare models we calculated four values for each model; Akaike’s Information Criterion (AIC), delta AIC, relative likelihood, and AIC weight. Smaller AIC values indicate “better” models and when comparing two models we calculated the difference in AIC values (delta AIC) (Akaike, 1973; Burnham & Anderson, 2002). A delta AIC of less than 2 indicates little difference between competing models; a delta AIC of 2–10 indicates moderate support for a difference between the models, and a delta AIC of greater than 10 indicates strong support (Burnham & Anderson, 2002). Relative likelihood represents the probability of a model given the data, whereas AIC weight is the discrete probability of each model (Burnham & Anderson, 2002).

The best model from each dataset was defined as having a delta AIC of 0.00, although preference was given to the model with the fewest variables if two or more models had a delta AIC of less than 2. This process was repeated for all six data sets: Puget Sound first, second, and third summer of growth and Strait of Georgia first, second, and third summer of growth.

Age, length, and growth

We estimated the age of 439 Pacific hake collected from PS and SoG in 1979–2016 (Table 1). These fish ranged in age from young-of-the-year to 6, with 2 year-olds being most abundant and comprising 38% of all individuals (Fig. 3A). Based on these ages and when fish were collected, we estimated hatch year to range from 1974 to 2015 (Table 2). Median total length (mm) at capture was 243 mm and 410 mm for PS and SoG, respectively, with SoG fish showing a bi-modal size distribution (Fig. 3B).

Our reconstruction of Pacific hake length at each summer and winter together with our repeated measures ANOVA, across all hatch years, indicated a significant interaction between population and year/season of growth (F = 339.7; df = 1,1950; p = 0.00022). By the second summer, individuals from SoG were, on average, 18% larger than fish from PS (mean size of 265 mm and 215 mm, and standard deviation of 42 mm and 28 mm, respectively) (Fig. 4A). The repeated measures ANOVA comparing estimates of length per year/season of SoG fish did not reveal significant differences between those that hatched in the 1970s versus 1990s (Fig. 4B). However, the same analysis performed on PS fish revealed a significant interaction between decade in which they hatched and year/season (F = 40.6; df = 3,1201; p = 0.00022). Tukey’s post hoc test showed significant differences between fish that hatched in the 2010s to those that hatched in both 1980s and 1990s (Fig. 4C).

Repeated measures ANOVA of our growth estimates, across all hatch years, indicated a significant interaction between population and year/season (F = 45.4; df = 1,1950; p = 0.00015). By the second summer individuals from SoG grew, on average, 42% more than fish from PS (mean growth 120 mm and 69 mm, and standard deviation of 28 mm and 31 mm, respectively) (Fig. 5A). The repeated measures ANOVA comparing estimates of growth per year/season of SoG fish did not reveal significant differences between those that hatched in 1970s versus 1990s (Fig. 5B). When we performed the same analysis on fish from PS it revealed a significant interaction (F = 15.7; df = 3,1201; p = 0.00031) between decade in which they hatched and year/season. Tukey’s post-hoc test showed significant differences among fish that hatched in the 1980s, 1990s, and 2010s, as well as between those that hatched in the 1980s and 2000s (Fig. 5C). The first and second summers showed the most dramatic increases and decreases, respectively, in growth across the decades in which fish hatched. Specifically, in their first summer fish that hatched in the 2010s grew, on average, 147 mm compared to 108 mm for those that hatched in the 1980s. However, fish that hatched in the 2010s grew 31 mm in their second summer while those that hatched in 1980s grew 107 mm (Fig. 5C).

Generalized linear model

We obtained elemental concentrations from 214 individuals collected in PS and 132 from SoG. All isotopes analyzed were found in concentrations above the detection limit for 100% of the otoliths sampled. We observed an acceptable level of analytical precision based on NIST (n = 87) mean coefficient of variation and standard deviation: Mg:Ca (mean = 2.35, standard deviation = 1.44), Mn:Ca (mean = 2.42, standard deviation = 1.09), Sr:Ca (mean = 1.65, standard deviation = 0.52), and Ba:Ca (mean = 2.29, standard deviation = 1.13).

We investigated variability in somatic growth separately for each population owing to the significant differences in growth and relatively little overlap in hatch years between populations. Overall, our GLM approach to understand what factors explain variability in somatic growth rate revealed a different suite of variables, more or less, with each summer and population. Further, AIC weights were relatively low (<0.3) indicating that a substantial amount of variability in growth was unexplained.

For SoG fish, results of the GLM analysis indicated five models that explained the greatest amount of variability in the first summer of somatic growth and were indistinguishable (i.e., delta AIC < 2.0) (Table 3). Of these five models the simplest model was chosen as the best model and it showed a negative relationship between somatic growth and Mn:Ca and Sr:Ca, and a positive relationship between somatic growth and Ba:Ca, ‘Harbor seal abundance’, and ‘Sea surface temperature’ (Figs. 6A–6E). Our GLM analysis of the second summer of growth revealed nine models that best explained variability and were indistinguishable (Table 3). Of these nine models, the one with the lowest AIC was also the simplest model. This model showed positive and negative relationships between somatic growth and Mn:Ca and Sr:Ca, respectively (Figs. 6F–6G). Finally, the GLM analysis of the third summer of growth showed that nine models best explained variability and were indistinguishable (Table 3). Of these nine models, three each consisted of only one variable: ‘Sea surface temperature’, ‘Harbor seal abundance’, and ‘Year’. Each of these three models showed a negative relationship to somatic growth (Figs. 6H–6J).

With respect to PS fish, the GLM analysis of the first summer of growth showed that three models best explained variability in somatic growth and were indistinguishable (Table 3). Of these, the best was also the simplest. This model indicated positive relationships between somatic growth and year, Mn:Ca, Ba:Ca, and ‘Harbor seal abundance’, and negative relationships with Sr:Ca (Figs. 7A–7E). Results of the GLM analysis of the second summer of growth revealed that nine models best explained variability in somatic growth and were indistinguishable (Table 3). Of these, the best was also the simplest, in which a positive relationship was observed between somatic growth and Mn:Ca and a negative relationship to Sr:Ca and ‘Harbor seal abundance’ (Figs. 7F–7H). Lastly, the GLM analysis of the third summer of growth showed that seven models best explained variability and were indistinguishable (Table 3). Of these, two were the simplest model, both consisted of four variables, and had Mn:Ca, Ba:Ca, and ‘Harbor seal abundance’ in common (Table 3). Both models showed a positive relationship between growth and Mn:Ca, and negative relationships with Ba:Ca and ‘Harbor seal abundance’ (Figs. 7I–7L). In one model there was a negative relationship between growth and year, and a lack of a relationship with PNI.

Size-at-age and growth differences between populations

Our findings of size-at-age differences between PS and SoG populations (Fig. 4A) were expected given that the listing of the Georgia Basin DPS as a Species of Concern was due, in part, to the reduced size-at-age of PS fish observed through the 1990s (Gustafson et al., 2000). For example, Pedersen (1985) and Goñi (1988) reported that PS fish collected in the mid-1960s to early 1980s grew approximately 2–5 cm shorter than those from SoG. What our study also revealed is that this divergence in body size between fish from PS and SoG occurred during their second summer of life (Fig. 4A). Size differences further increased by their third summer when PS fish grew, on average, 73% less than fish from SoG (i.e., 23.2 mm compared to 88.7 mm, respectively; Fig. 5A). It should be noted, however, the limited temporal overlap of fish that hatched in PS and SoG. Specifically, only the hatch years from 1994 through 1998 included individuals that were sampled from both PS and SoG (Table 2).

From 1965–1974 most of the PS fish ranged in size from 32 cm to 45 cm (average length of 36.2 cm, with recorded maximum lengths of 45 cm and 73 cm for males and females, respectively) (Pedersen, 1985), yet since 2008, less than 5% of sampled PS fish (n = 11,652) had a length greater than 30 cm (Blaine, Lowry & Pacunski, 2020) and Supplemental Table 1). Pedersen (1985) suggested that the declining body size in PS through time, and the fact that PS fish mature at a smaller size, relative to SoG fish, may have been due to the intense commercial fishery in PS in the 1970s and 1980s. Specifically, the extensive removal of large individuals from a population has been shown to select for individuals that mature earlier and grow less (Allendorf & Hard, 2009; Olsen et al., 2004; Sharpe & Hendry, 2009). For instance, California sheephead (Semicossyphus pulcher) had significantly smaller body size, earlier maturation, and a reduced maximum lifespan following more than 20 years of high fishing pressure (Hamilton et al., 2007). Similarly, zebrafish (Danio rerio) invested more in reproduction and attained a smaller adult body size when exposed to intense size-selective harvesting (Uusi-Heikkilä et al., 2015). Interestingly, models simulating zebrafish population dynamics indicated a slow population recovery even after a moratorium on harvest. Uusi-Heikkilä et al. (2015) suggested that this lack of population recovery was linked to changes in body size and reproductive investment, which placed the population at a disadvantage to manage natural selection pressures that often favor large body size. Similar harvest-related changes to life-history could explain why the smaller body size and the lack of recovery of the PS population persists today despite the closure of the fishery for thirty years.

Size-at-age and growth differences through time and within populations

We observed significant differences in size-at-age of PS fish such that a shift was seen towards smaller fish through time (Fig. 4C), yet no differences in size-at-age were detected between SoG fish that hatched in the 1970s and 1990s (Fig. 4B). In terms of somatic growth, we found a similar pattern in both populations whereby in recent years fish tended to grow more in their first summer and less in subsequent summers. Specifically, SoG fish that hatched in the 1990s grew, on average, 11% more in the first summer, and 25% less in the third summer, relative to those that hatched in the 1970s (Fig. 5B). In comparison, PS fish that hatched in 2010s grew 26% more in their first summer than those that hatched in 1980s, yet 71% less in the second summer compared to those that hatched in 1980s (Fig. 5C). Beamish & McFarlane (1999) suggest that this increased somatic growth during the first summer is linked to a shift in climate and ocean conditions that occurred in the late 1980′s. This regime shift is hypothesized to have caused a greater temporal overlap between young-of-year hake and an important prey item (e.g., copepods) in SoG (King & McFarlane, 2006). Because PS and SoG populations are part of the same basin and share similar characteristics it is possible that the same processes played a role in the size and growth patterns in SoG would be seen in PS.

Our observed decrease in growth during the third summer of SoG fish is consistent with findings by King & McFarlane (2006) whereby an increase in hake abundance following the aforementioned regime shift likely lead to reduced somatic growth later in life because of increased density dependent competition. The third year of growth is an important period in hake development in which they become reproductively mature and are suspected to move deeper to feed (Saunders & McFarlane, 1999). Therefore, if intraspecific competition is high during the time when surplus energy is allocated to reproductive development then it is expected that reduced energy will be available for growth. Regardless of the mechanism(s), the declines in growth through time, in both populations, highlights increased concern for the sustainability of this Species of Concern.

Our investigation of what factors and/or processes explain variability in seasonal growth indicated some similarities between populations. For instance, growth in the first summer was best explained by models that included positive relationships to harbor seal abundance and to barium levels (i.e., Ba:Ca) (Table 4). The positive growth associated with harbor seal abundance could indicate that hake resources are density-mediated and therefore predation-induced reductions in hake abundance may have led to reductions in intraspecific competition (Heithaus et al., 2008). Ba:Ca is positively associated with areas of nutrient rich upwelling (Mohan et al., 2018; Wheeler et al., 2016) and freshwater input (Elsdon & Gillanders, 2005; Hamer, Jenkins & Coutin, 2006; Miller, Gray & Merz, 2010), and the positive relationship between Ba:Ca and growth (Figs. 6C and 7D) might be due to increases in productivity, which could have increased prey (e.g., copepods) abundance due to nutrient enriched waters brought to the surface and/or from riverine sources. For instance, Takahashi et al. (2012) reported a positive relationship between otolith growth of northern anchovy (Engraulis mordax) and upwelling intensity, which they attributed to increased productivity and in turn improved prey nutritional value. Strontium (i.e., Sr:Ca) in otoliths is reported to be positively associated with salinity (Macdonald & Crook, 2010; Martin & Thorrold, 2005; Zimmerman, 2005). Thus the negative relationship we observed in both populations between growth (in the first and second summers) and Sr:Ca (Table 4; Figs. 6B, 6G, 7C, and 7G), and the aforementioned relationship to Ba:Ca, suggests that the first two years of growth are sensitive to freshwater flow, such that growth declines with a reduction in freshwater. Main spawning aggregations for both PS and SoG populations occur in close proximity to major sources of freshwater (e.g., Stillaguamish and Snohomish rivers in PS and Fraser River in SoG). This freshwater input would cause a stratified layer of low-salinity water above the well-mixed marine layer. This is significant since Pacific hake larvae are known to aggregate in surface waters (Gustafson et al., 2000).

The third summer of growth for fish from both populations was negatively associated with harbor seal abundance, which is an opposite pattern observed for the first summer of growth. These positive and negative relationships may suggest that younger hake experience a competitive release in years with more harbor seals, while older hake experience reduced performance possibly due to energy expenditure through predatory avoidance. Heithaus et al. (2008) noted that aside from direct mortality, predators can strongly influence prey behavior such as foraging. It is likely that hake recovery, particularly in PS, is limited given that harbor seal abundances increased by 7% (Saunders & McFarlane, 1999) and 3% (Gustafson et al., 2000) in SoG and PS, respectively, in the early to mid-1990s and appear to be stable in both populations (Jeffries et al., 2003; Olesiuk, 1999).