A Dynamic Stress-Scape Framework to Evaluate Potential Effects of Multiple Environmental Stressors on Gulf of Alaska Juvenile Pacific Cod

Introduction

Environmental variability and associated trends induced by climate change can create unique and potentially stressful conditions for marine fisheries that fluctuate through space and time. Increasing water temperatures and ocean acidification related to changing ocean conditions have the potential to negatively impact individual fish – through direct and indirect effects – which can potentially impact fisheries by way of fishable biomass, annual stock productivity, and spatial shifts in the population outside of the traditional fishing areas (Holsman et al., 2020; Laurel and Rogers, 2020). The Gulf of Alaska Pacific Cod (Gadus macrocephalus) fishery is one example of a commercial fishery that is highly susceptible to changing ocean conditions. In recent years, Pacific Cod abundance has declined dramatically in the Gulf of Alaska (GOA) ecosystem, leading to a fisheries disaster declaration in 2018 and a closure of the federal fishery in 2020 (Barbeaux et al., 2020). Quantifying the spatial and temporal footprint of modern environmental stressors on GOA Pacific Cod is crucial for ensuring the accuracy of plans and forecasts for this valuable fishery.

Decadal and long-term patterns influence the GOA’s thermal variability due to warm-phase shifts of climatic phenomena and background warming related to climate change. The GOA has also experienced two anomalous marine heatwave events in the past decade. Between late 2013 and 2016, a marine heatwave occurred in the GOA that exceeded the magnitude and duration of any other heatwave on record in the region. This heatwave event led to temperature anomalies greater than 2.5°C (Bond et al., 2015; Di Lorenzo and Mantua, 2016) and unprecedented shifts in the region’s biological communities, including increases in harmful algal blooms, reductions in fishery recruitment, and mass mortality of marine mammals and seabirds (Leising et al., 2015; Piatt et al., 2020; Suryan et al., 2021). In the summer of 2019, the GOA was impacted by another marine heatwave, leading to similarly extreme temperature anomalies (Amaya et al., 2020). The increasing frequency and intensity of marine heatwaves in the coming years are likely to coincide with challenges to marine fisheries’ long-term health (Oliver et al., 2018; Cornwall, 2019).

In addition to thermal variability, there are several sources of environmental variability that co-occur in the GOA and which have the potential to overlap, compound, or negate one another in ways that are still not fully understood (Bopp et al., 2013). Ocean acidification is one such potential stressor and refers to the gradual increase in hydrogen ions and consumption of carbonate ions in marine systems in response to elevated levels of anthropogenic carbon dioxide entering the ecosystem (e.g., Feely et al., 2004; Evans et al., 2014). The GOA, with the naturally lower concentrations of carbonate ions characteristic of high latitude marine habitats, is at higher risk for the effects of ocean acidification as additional losses of carbonate ions can lead to proportionally larger changes in seawater chemistry than would be expected in lower latitude systems (Fabry et al., 2009; Mathis et al., 2011). Additionally, the GOA and Kodiak Island regions are known to have higher seasonal variability in natural surface pCO₂ levels (seasonal amplitudes of as much as 309.8 and 279.4 μatm, respectively) when compared to open ocean environments due to upwelling conditions (Chen and Hu, 2019). Some interactions between sea surface temperature and the carbonate system are known – for example, carbon dioxide becomes less soluble in water as water temperature increases – but the full effects are not yet known in the context of physiological stress to living marine resources (Mathis et al., 2015).

Recent declines in Pacific Cod stocks are likely connected to changing thermal conditions in the GOA. Pacific Cod is a stenothermic species with peak hatch success occurring in 4–5.5°C water (Laurel and Rogers, 2020) and optimal juvenile growth occurring below 13°C in nursery habitats (Hurst et al., 2010). It is likely that Pacific Cod, like other stenothermic species, exhibit stage-specific responses to warming, with potential for thermal bottlenecks between different life stages (Dahlke et al., 2020; Rogers et al., 2020). Early life stages may be particularly sensitive to changes in temperature due to higher metabolic rates in the embryonic, larval, and juvenile stages (e.g., Finn et al., 2002: Werner, 2002). Elevated temperatures in the GOA are associated with high larval mortality of Pacific Cod (Doyle and Mier, 2016) and a larger size at nursery entry. During the 2014–2016 and the 2019 marine heatwaves, it is hypothesized that a loss of suitable spawning habitat during heatwave conditions contributed to declines in Pacific Cod hatch success and consequent older life stages (Laurel and Rogers, 2020).

Pacific Cod also appear to be moderately impacted by ocean acidification in the GOA. Elevated levels of carbon dioxide in the water column during the first two weeks of the larval duration are associated with suppressed growth, although this trend has been shown to reverse as larvae grow (Hurst et al., 2019). By the juvenile stage, effects of elevated carbon dioxide levels on growth largely disappear (Hurst et al., 2012). It is hypothesized that fishes are better able to withstand the effects of ocean acidification than calcifying marine invertebrates due to their high metabolism and ability to regulate their inter/intracellular acid-base balance (e.g., Melzner et al., 2009). However, elevated levels of carbon dioxide could lead to several indirect impacts on fishes including behavioral changes (e.g., Williams et al., 2019); changes in the biochemical structure of prey (e.g., Cripps et al., 2016); and changes in the abundance of calcifying marine invertebrate prey (e.g., Lischka et al., 2011). Many of the direct and indirect impacts of ocean acidification on cod and other marine finfish species are unknown, and the modest impacts of ocean acidification observed in the laboratory may indeed be “significant” when compounded by other environmental stressors or realized over longer time scales. Here, we do not account for potential interaction between temperature and ocean acidification as further studies are needed to better understand the interactive effects of these two variables on Pacific Cod early life stages.

Quantifying and visualizing environmental conditions in space and time can be beneficial in understanding when and where the Pacific Cod fishery may be most influenced by changing ocean conditions in the GOA. By viewing the ocean as a dynamic mosaic of physical and biogeochemical seascapes, researchers, managers, and other end-users can efficiently communicate information about environmental conditions in the ocean through space and time. This approach is known as dynamic seascape classification (Kavanaugh et al., 2014). Recent advances have allowed for the dynamic classification of ocean regions based on multiple biophysical interactions observable by satellite remote sensing or marine ecosystem models (e.g., Kavanaugh et al., 2018). These dynamic classifications have been used to predict biogeochemical processes (Kavanaugh et al., 2014), phytoplankton structure and function (Kavanaugh et al., 2015; Montes et al., 2020), and individual species occurrences (Breece et al., 2016). The seascape ecology framework has potential utility for quantifying biological and human responses to known environmental stressors, including elevated temperature and ocean pCO2. By applying a dynamic seascape framework in marine resource management contexts, Klajbor (2020) established possible novel relationships between dynamic ocean management principles and socioeconomic vulnerability in fishing communities.

Large, biophysical time series like those used in a seascape ecology approach have been classified in space and time simultaneously using a probabilistic self-organizing map (prSOM) machine learning algorithm (Kavanaugh et al., 2014, 2016, 2018). The prSOM combines two statistical algorithms, K-Means and expectation-maximization, and a machine learning algorithm called the Self-Organizing Map (Kohonen, 1990; Richardson et al., 2003) to approximate the multivariate distribution of data as a mixture of Gaussians and uniform distributions (Anouar et al., 1998; Lebbah et al., 2015). The prSOM algorithm has utility in capturing biological metrics, such as Pacific Cod growth, in addition to environmental metrics because it is robust to non-linear relationships between variables and is topology-preserving. Also, the prSOM algorithm does not rely on data being identically distributed or independent (i.i.d). This feature is critical for classifying complex, potentially non-linear interactions like those found in data from spatio-temporal time series, which are not independent because the current year’s classification almost always depends on previous months, years, and nearby regions.

Here, we introduce “stress-scapes” as a novel machine learning framework to understand juvenile Pacific Cod’s response to potential environmental stressors across the GOA ecosystem. Juvenile Pacific Cod was selected for this project due to the availability of information and data related to the effects of temperature and ocean acidification on early life stages. While there remain many unknown biological responses to these stressors, both alone and in their potential interactive effects, we incorporate the most recent available data on this critical species for this proof-of-concept study. Stress-scapes incorporate remotely sensed sea surface temperature (SST) data, pCO₂ data, and estimated juvenile Pacific Cod specific growth rate (SGR) into a prSOM algorithm to identify ocean conditions that may be stressful to juvenile Pacific Cod across a broad spatial and temporal scale. The juvenile Pacific Cod stress-scape presented in this study is not intended as a comprehensive reconstruction of Pacific Cod growth history in response to environmental variability in the GOA, but rather as a novel framework to visualize and classify potential environmental stress on the species across space and time. In this way, the stress-scape framework is distinct from an Individual-Based Model (IBM), which is a biophysical model that describes fish growth, recruitment, and mortality as dependent on environmental factors (Hinckley et al., 1996; Hermann and Moore, 2009; Koenigstein et al., 2016). IBMs are often coupled with Regional Ocean Modeling Systems (ROMS) and dispersal models to capture environmental conditions at depth and track the movement of fish in early life stages through their environment (e.g., Stockhausen et al., 2019). An IBM for Pacific Cod has already been developed for the GOA ecosystem (Hinckley et al., 2019).

In our case study, we develop weekly stress-scapes from 2010 to 2016, including marine heatwave conditions in 2014–2016, using existing juvenile Pacific Cod laboratory growth models and remotely sensed environmental data. These stress-scapes are then evaluated spatiotemporally to determine how well they capture changing conditions both seasonally and interannually, including during the marine heatwave event. We then evaluate the sensitivity of the stress-scape model to the Pacific Cod SGR input variable to determine the importance of this variable in influencing the model’s classification. Finally, we compare the stress-scape model output to otolith-derived growth data from juvenile Pacific Cod collected from Kodiak Island, AK to comment on the success of the model in capturing processes occurring in the field.

Materials and Methods

Environmental Data

Sea surface temperature and sea-surface carbonate chemistry data (pCO₂) were collated for the GOA ecosystem between the years 2010–2016 (see Supplementary Material: Appendix A for spatial extent of environmental data). Environmental data were selected from publicly available modeled, in situ, and remotely sensed products. The data’s temporal range was selected to include dates before, during, and after the 2014–2016 marine heatwave event. Sea surface temperature data were subset from the National Oceanic and Atmospheric Administration (NOAA) Optimum Reynolds Interpolation SST V2 product provided by the OAR/ESRL PSL, Boulder, Colorado, United States (Reynolds et al., 2002). The product is developed using a combination of in situ and remotely sensed data, producing weekly means on a one-degree grid. Associated variance related to each point was also downloaded and used to calculate total error related to the stress-scape classification. Sea-surface carbonate chemistry data were subset from the NOAA National Center for Environmental Information CSIR-ML6 v2019a Global Sea-Surface pCO2 product. A smoothed mean monthly pCO₂ product was selected from their ensemble mean product that is calculated using a combination of 6 machine learning techniques. The smoothed product differs from their raw product, as a 3 × 3 × 3 convolution is applied to the data through time, latitude, and longitude (Gregor et al., 2019). Mean monthly values ranged from 262.6 to 469.4 μatm in the study domain.

Only the environmental data from the known spatial and temporal ranges of Pacific Cod were included in the subsequent machine learning training and stress-scape visualization, ensuring that the selected data and resulting stress-scapes were relevant to juvenile Pacific Cod distribution. Age-0 juvenile Pacific Cod are distributed in shallow, nearshore habitats, typically <30 m depth during their first summer (Laurel et al., 2007, 2009). We used a bathymetric data subset from the General Bathymetric Chart of the Oceans 2020 Gridded Global Bathymetry product (GEBCO Bathymetric Compilation Group 2020, 2020) for 30 m and 100 m isobaths to compare habitat size across both depth distributions. Due to the small size of the habitat at <30 m depth, we selected 100 m depth as the cutoff for the juvenile habitat, acknowledging that this represents a broader habitat range than would be expected for a typical age-0 juvenile Pacific Cod during their first summer and subsequently may bias our interpretation toward cooler temperatures than would be experienced in situ (e.g., Laurel et al., 2012). We used a year-round temporal range for the juveniles in the nursery habitat to represent both the first summer in the nursery and overwintering in the nursery. The study’s spatial domain only includes the GOA, so any data that our filter included from the Bering Sea were excluded from the training of the prSOM model. Additionally, there were 5 points where the pCO₂ or SST were not available (Supplementary Material: Appendix B); rather than attempt to interpolate the data where it was not possible due to missing data, the value of each pixel was set to NaN.

To account for seasonal variance the monthly z-scores were computed for each grid point across the GOA for each of the relevant variables:

where j is the variable index, i is the observation index and the represents an average value across the entire GOA through time. The spatio-temporal z-scores were used to train the prSOM to compute stress-scape classes across the GOA from 2010 to 2016 (Gregor et al., 2019); hence, the stress-scape classes represent regions that similarly deviate from the spatio-temporal mean.

Juvenile Pacific Cod Growth Data

Juvenile Pacific Cod temperature-dependent growth data were obtained from two existing experimental studies: Laurel et al. (2016) and Hurst et al. (2010). These studies were all conducted in the laboratory using age-0 juveniles collected from the GOA population. Similar growth data related to variable pCO₂ are currently unavailable for age-0 juvenile Pacific Cod [although see Hurst et al. (2019) for earlier life stage responses]. As a proxy, we used the pCO₂ growth response of age-0 juvenile Walleye Pollock (Gadus chalcogrammus) based on experimental work from Hurst et al. (2012). While response to ocean acidification is variable among marine finfish, Walleye Pollock are congeners of Pacific Cod and share a similar growth response at the age-0 juvenile stage (Laurel et al., 2016).

A weight-based growth model was used instead of a length-based growth model throughout this stress-scape analysis to better align with potential future bioenergetic applications of the stress-scape to fish bioenergetics. Therefore, a specific growth rate (SGR) model was used to approximate the percent change of fish wet weight over time. SGR were calculated from the formula:

where g is the instantaneous growth coefficient obtained from the equation:

where WW₁ is the wet weight of the fish at the beginning of the experiment at time t₁, and WW₂ is the wet weight of the fish at the conclusion of the experiment at time t₂.

A temperature-dependent growth model was developed using a general linear model with SGR as the independent variable and a quadratic temperature term as the dependent variable (Simple Linear Regression; R²_adj = 0.636; P < 0.0001). Uncertainty in the model was evaluated using 95% confidence intervals and 95% prediction intervals.

Juvenile Pacific Cod growth response to ocean acidification was divided into two potential outcomes: growth under “low” pCO₂ conditions (<1,000 μatm) and growth under “high” pCO₂ conditions (>1,000 μatm). These two classifications represent growth at ambient to moderate pCO₂ concentrations (“low”), and growth at high to very high pCO₂ concentrations (“high”), because juvenile Pacific Cod growth only begins to change as pCO₂ concentrations become very high (e.g., Walleye Pollock; Hurst et al., 2012). Growth under low pCO₂ conditions was unchanged from the temperature-dependent growth curve and was represented by the equation:

Growth curves under high pCO₂ conditions were calculated from the percent decrease in growth (−2.3%) observed in Hurst et al. (2012). Specific growth rates from the temperature-dependent growth curves were then depressed to reflect the percent change in growth that would be expected under elevated pCO₂ conditions. Growth under high pCO₂ conditions was represented by the equation:

Classification of Stress-Scapes

Stress-scapes were classified using remotely sensed observations of SST; partial pressure of CO₂ (pCO₂) as an indicator of the carbonate system; and the specific growth rate (SGR) of juvenile Pacific Cod. The prSOM approximates the multivariate distribution of SST, pCO₂ and SGR through space and time. The prSOM was configured with 225 neurons in a 15 × 15 neuron map. Previous efforts in the North Pacific found that this configuration resulted in sufficiently large number of data points (e.g., >500) per node (Kavanaugh et al., 2014). The prSOM algorithm executed 50 epochs in order to give the prSOM enough time to converge. We applied Hierarchical Agglomerative Clustering (HAC) to the 225 neurons grouping them into five individual stress-scapes (e.g., Saraceno et al., 2006; Hales et al., 2012; Kavanaugh et al., 2018).

The SST, pCO₂, and SGR maps of the GOA were interpolated using bilinear interpolation from 14 × 40 grid points to a 780 × 2340 image to classify stress-scape class boundaries in the remotely sensed data. The interpolated image was copied to the Graphics Processing Unit (GPU) memory and classified according to the trained neurons’ weight vectors and variances by identifying the stress-scape class most likely to generate each interpolated value in the GOA. Our unique prSOM implementation classified the GOA in parallel using the GPU. Parallel classification is an important distinction because many visualization methods rely on processing done on the Central Processing Unit (CPU), but the amount of data we needed to classify made classification on the CPU impractical. Any consumer-grade GPU released in the past decade can classify the entire GOA in a small number of clock cycles (∼1–100 clock cycles) by taking advantage of massively parallel processing. In contrast, even with a parallel CPU implementation, the number of clock cycles will range from 300,000 to 1.7 million clock cycles on most modern processors. There are 365 time slices, so the classification must be done 365 times; hence the speedup is substantial.

We selected 5 stress-scape classes from the HAC using the Gap Statistic (Tibshirani et al., 2001) which suggested 5 classes with 2 well separated classes and 3 less separated classes. Choosing an optimal number of stress-scapes is difficult as there is no single metric that can select an ideal number of clusters with certainty, and often, the statistically significant number of clusters does not agree with domain expert analysis. For our implementation, we compared results from three different methods of clustering evaluation: the elbow method, the silhouette method, and the Gap Statistic (see Supplementary Material: Appendix C). The elbow method works by estimating the amount of variance in the data explained by the clusters, an optimal number of clusters is chosen by identifying the point in the graph where an increase in the number of clusters does not provide additional explanation of the variance. The silhouette method works by estimating the quality of each cluster and then reporting the average quality over all the clusters for increasing number of clusters. An optimal number of clusters can be determined from the silhouette method by choosing a number K for which the quality of the clusters is greater than the quality of the K-1 and K + 1 clusters. The Gap Statistic works by computing the HAC algorithm for an increasing number of clusters. For each number of clusters, from 1 to 10, we computed the inter-cluster sum of square error, the average quality of the clusters, and the gap statistic using the clustGap() function in R.

In order to evaluate the effects of anomalous heatwave conditions on the stress-scape classification results, classification results from 2010 to 2013 were compared to classifications during a known anomaly, the marine heatwave event of 2014–2016, to more normal ocean conditions preceding the heatwave. In comparing the classification results, we identified a dominant stress-scape and examined the qualities of the dominant stress-scape using a decision tree and a table to records the frequency and the environmental characteristics of the stress-scapes. Further details describing our implementation of the prSOM are provided in Supplementary Material: Appendix C.

Sensitivity Analysis

We trained a second, bivariate prSOM to compare the results of the stress-scape to a classification that did not include juvenile Pacific Cod SGR. This bivariate classification allowed us to evaluate the sensitivity of the trivariate stress-scape framework to the additional information gained by including the juvenile Pacific Cod SGR variable in the training process. To account for the reduced variance in the bivariate prSOM, we used the Gap Statistic to estimate 4 classes in the bivariate prSOM. We applied Hierarchical Agglomerative clustering until the 225 neurons were merged into 4 distinct stress-scape classes. We focused our comparison efforts on the coastal regions near Kodiak Island.

Pacific Cod Field Growth

Otolith structural analysis of juvenile Pacific Cod collected from Kodiak Island nurseries is being completed as part of a related study (Thalmann et al., unpublished data). Therefore, after developing the lab-based model for juvenile Pacific Cod SGR, training the prSOM algorithm, and generating the mosaic of stress-scapes across the coastal GOA, we estimated otolith-based SGR from these field-collected juveniles to compare with our stress-scape classifications of growth observed in nearshore Kodiak Island nursery habitat. This comparison aimed not to provide a comprehensive analysis of otolith data but rather to identify a relationship between temporal patterns in the otolith-based growth rates and temporal patterns in the stress-scapes. Juvenile Pacific Cod were collected in annual summer beach seine surveys on Kodiak Island in 2010 and 2012–2016 through the NOAA Alaska Fisheries Science Center Fisheries Behavioral Ecology program. Samples were not available in 2011. Pacific Cod were collected from two shallow (<30 m) embayments along the northeastern coast of Kodiak Island during successive days in July each year. Fish were captured using a demersal beach seine (36 m bag with 5-mm mesh; 1 m × 2.25 m wings with 13-mm mesh), sorted to species, counted, and measured (total length, mm). Samples were then frozen (−20°C) and returned to Hatfield Marine Science Center in Newport, OR for post-processing. During post-processing, standard length (mm) and wet weight (g) were measured, and otoliths were removed for analysis (Table 1).

Table 1. Annual sample sizes, average fish size, mean July ocean temperature, mean July pCO2 levels, and ocean conditions for juvenile Pacific Cod captured annually in July 2010 and 2012–2016 from Kodiak Island, AK.

Otolith structural analysis was used to generate growth estimates for field-collected juveniles collected at Kodiak Island. Otoliths are calcium carbonate structures in the inner ear of teleost fishes that can be used to determine age and growth of an individual (Pannella, 1971; Campana and Thorrold, 2001). Otolith size and body size are highly correlated in Pacific Cod (DiMaria et al., 2010; Miller et al., 2016), and daily formation of otolith increments has been validated for Pacific Cod at 10°C up to 120 days post-hatch (Narimatsij et al., 2007). Sagittal otoliths were mounted on glass slides using thermoplastic resin and polished to expose the core in the transverse plane using Wetordry paper (800–2000 grit), Buehler lapping film (3–30 micron grit), and alumina slurry (0.3 micron). Polished otoliths were imaged at 100× and 400× magnification using a Leica DM1000 compound microscope and a Levenhuk M1000 digital camera. Final otolith radius and otolith radius from 8 and 16 days prior to capture were quantified using increment analysis. Daily increments for these two 8-day growth periods (representing three weeks of nursery residence at the same 8-day resolution as the stress-scape model) were counted and measured using ImagePro Premier software. Each otolith was interpreted three times with greater than 80% similarity in mean increment width between reads, and reads were averaged in order to compare otolith-generated estimates of growth to the stress-scape.

First, we estimated fish length and mass at 8-d and 16-d prior to capture based on otolith increment analysis, and then estimated growth in length and mass over these two 8-d periods for comparison with the 8-d prSOM output. Otolith radius measurements and fish standard length were log-transformed to meet parametric assumptions, and a log-linear model with otolith radius, year, and their interaction as independent variables and standard length as a dependent variable was developed (Multiple Linear Regression; R²_adj. = 0.85; P < 0.0001; Supplementary Material: Appendix D). Using this relationship, we developed a proportional back-calculation model to estimate fish size 8- and 16-d prior to capture using the measured otolith radius from 8- and 16-d prior to capture (Francis, 1990; Campana and Jones, 1992). Length-based growth of juvenile Pacific Cod in mm per day was examined for each year and compared to a mm/d laboratory growth model obtained from the same studies (Hurst et al., 2010; Laurel et al., 2016) as the laboratory SGR model used in the stress-scape analysis.

To estimate juvenile wet weight 8- and 16-d prior to capture, we developed a second log-linear model with standard length and year as independent variables and wet weight as a dependent variable (Multiple Linear Regression; R²_adj. = 0.96; P < 0.0001; Supplementary Material: Appendix D). Wet weight from 8- and 16-d prior to capture was then determined using a second proportional back-calculation model using estimated standard length from 8- and 16-d prior to capture. Juvenile Pacific Cod SGR in the field was calculated in the same manner as the laboratory SGR model described above from Hurst et al. (2010) and Laurel et al. (2016). We utilized weight-based SGR calculated from otoliths rather than length-based growth estimates for comparison to the stress-scape.

To compare otolith-generated SGR estimates from the 8- and 16- day periods prior to capture with the stress-scape, we determined the dominant stress-scape class for coastal Kodiak Island at an 8-day resolution in July for each year between 2010 and 2016. Mean SGR from the dominant stress-scape class was then compared to otolith-generated SGR for that year using a paired t-test. Otolith-generated growth estimates for each year were additionally compared to the temperature-dependent laboratory growth model using a paired t-test. We also compared otolith-derived growth rates in heatwave conditions in 2014–2016 to more normal conditions in 2010–2013 using a 1-way analysis of variance.

Results

Stress-Scape Analysis

Utilizing the prSOM and HAC, we resolved 5 unique stress-scapes across the GOA that quantified the spatio-temporal footprint of changing ocean conditions and the modeled specific growth rate (SGR) of juvenile Pacific Cod. Table 2 describes the means and standard deviations for the SST, pCO₂, and SGR values of each of the resolved classes. Stress-scape classes were determined and distinguished semantically through a decision-tree of z-scores for each input variable, including the characteristics of each decision mode, which effectively describes how the classes were sorted (Figure 1). Stress-scape classes 2, 3, 4, and 5 (maroon, green, orange, and yellow) were distinguished by their respective growth rate input variables when they were decided by node 5. This highlights the non-linear relationship between SST and SGR. Class 1 (blue) was primarily distinguished by differences in SST and pCO₂.

Table 2. Mean values for each of the 5 resolved stress-scapes.

Figure 1. The results of the decision tree trained on the output of stress-scape classification using 5-stress-scapes and 4 variables: SST z-score (ZSST), juvenile Pacific Cod SGR z-score (ZGROWTH) measured in percent wet weight per day (%WW/day), pCO₂ concentration (μatm) z-score (ZPCO2), and the month of the remotely sensed observation (MONTH). Each node of the decision tree represents a decision about the variable labeling the vertex and is labeled with a number. The root node represents a decision about the monthly z-score of SST. If the pixel’s SST is greater than.562 standard deviations above the mean SST value across the entire GOA from 2010 to 2016, the left branch is taken; otherwise, the right branch is taken. Note that by including SGR (ZGROWTH) as a classification variable we can make biologically relevant interpretations of the stress-scape classes that are classified by the right branch of node 1. The bar graph provides the percent of each class that is classified by each leaf node of the tree.

The stress-scape analysis identified distinct classes that corresponded to heatwave conditions in 2014–2016 compared to more normal ocean conditions in 2010–2013 for both coastal Kodiak Island (Figure 2) and the larger GOA coastal ecosystem (Figure 3). Animated versions of these figures with corresponding ribbon plots that show how the stress-scapes shift and progress through time can be found in Supplementary Material: Appendix E, Supplementary Videos 1, 2. At both the Kodiak and regional scales, seasonal and interannual patterns emerge in the incidence of each of the classes. Classes 2 and 5 contract over time, supplanted by expansions (in space and time) of classes 1 and 4. Class 2 is the low SST class that seems to be replaced by high pCO₂ in class 4. Class 5, the pre-heatwave late spring-early summer class, is replaced by class 3, characterized by high SST, during heatwave years. Class 3 starts appearing earlier and persisting each year during the heatwave.

Figure 2. Dominant yearly stress-scape classes for coastal Kodiak Island during the first week of July 2010 and 2012–2016. Heatwave conditions are present in 2014–2016. Column 1 highlights the dominant stress-scape classes for coastal Kodiak Island. Column 2 is the z-score for SST; Column 3 is the z-score for SGR, and Column 4 is the z-score for pCO₂. The red line marks the time shown in the figure across all years. Dark purple z-scores are below average and yellow z-scores are above average for that variable. An animated version of this figure can be found in Supplementary Material: Appendix E and Supplementary Video 1.

Figure 3. Ribbon plot showing total Stress-scape class area per time slice. Each column represents a single slice through time of the entire GOA. A stress-scape dominates in time if it persists along the x-axis. For example, in 2016, stress-scape class 1 was spatiotemporally dominant. From 2014 to 2016 stress-scape classes 1 and 4 dominate for a longer period than in 2010–2013. Classes 1 and 4 are coincident with the anomalous elevated growth rates seen in the otoliths data. A seasonal and interannual pattern emerges in the incidence of each of the classes: The spatiotemporal footprint of classes 2 and 5 contract over time, supplanted and are replaced by expansions of classes 1 and 4. Class 2, a low SST class, that seems to be replaced by high pCO₂ class 4. Class 5, the pre-heatwave summer class, is replaced by the high SST class 3, which starts appearing earlier and persisting each year during the heatwave. This ribbon plot can be viewed as a video alongside a map of the GOA in Supplementary Material: Appendix E and Supplementary Videos 1, 2.

In July, stress-scape class 5 (yellow) dominated the areas around Kodiak Island in 2010 and 2012, whereas stress-scape class 3 (green) expands both around Kodiak and within the greater GOA area during heatwave conditions in 2014, 2015, and 2016. In 2013, stress-scape class 3 and stress-scape class 5 were both present throughout the coastal GOA. Stress-scape class 3 was characterized by above average temperatures (10.42 ± 2.48°C), higher pCO₂, concentrations (359.4 ± 33.66 μatm), and high juvenile Pacific Cod SGR (1.40 ± 0.16% Wet weight/day). In contrast, stress-scape class 5 was characterized by near average temperatures (7.85 ± 2.25°C), lower pCO₂ concentrations (339.1 ± 36.41 μatm), and near average SGR (1.22 ± 0.23% Wet weight/day).

Sensitivity Analysis

The bivariate prSOM, trained on SST and pCO₂ but excluding SGR, resolved four stress-scape classes for the coastal GOA. Table 3 describes the means and standard deviations for the SST and pCO₂ values of each of the resolved classes used for the sensitivity analysis. By excluding the SGR variable, the dominant stress-scape for coastal Kodiak Island was reclassified to a new stress-scape class 1 (blue), which identifies a mean SST of 10.42 ± 2.48°C and mean pCO2 concentration of 359.4 ± 33.66 μatm across all years of analysis. Stress-scape class 1 of the bivariate model corresponds to stress-scape class 3 of the trivariate-model. Figure 4 shows the dominant yearly classes for the bivariate prSOM sensitivity analysis near coastal Kodiak Island during the first week of July 2010 and 2012–2016 and should be evaluated using Table 3. Stress-scape class 1 is distinguished from the other three classes with its higher-than-average SST and near average pCO₂ concentrations. The bivariate prSOM groups the two higher temperature classes of the trivariate dynamic stress-scape together because it lacks the SGR variable necessary for distinguishing between these classes. In addition, the bivariate prSOM is unable to distinguish between heatwave conditions and more normal ocean conditions in the GOA, with stress-scape class 1 dominant across all six years of the analysis.

Table 3. Mean values for each of the 4 resolved stress-scapes that were classified to test the significance of the growth rate variable.

Figure 4. Dominant yearly classes for the bivariate prSOM sensitivity analysis near coastal Kodiak Island during the first week of July 2010 and 2012–2016. Normal ocean conditions are present in 2010–2013, and heatwave conditions are present in 2014–2016. Column 1 highlights the dominant stress-scape classes for coastal Kodiak Island. Column 2 is the z-score for SST, and Column 3 is the z-score for pCO₂. The final column gives the areal extent of each stress-scape class. The red line indicates the time slice of the areal extent (see Supplementary Material: Appendix E and Supplementary Videos 1, 2).

Juvenile Nursery Growth

The dynamic stress-scape underestimated the observed juvenile Pacific Cod SGR in its classification schema for coastal Kodiak Island by 39–187% across all six years of the analysis (Figure 5). Mean July SGR from the dominant stress-scape for coastal Kodiak Island differed significantly from otolith-derived SGR for both the last 8 days of growth and the penultimate 8 days of growth (paired t-test; last 8 days: t₅ = 6.72; P = 0.0011; penultimate 8 days: t₅ = 6.94; P = 0.0009). Otolith-derived SGR were between 10 and 105% higher than the laboratory growth model across all six years for both the last 8 days of growth and the penultimate 8 days of growth (paired t-test; last 8 days: t₅ = 4.69; P = 0.0053; penultimate 8 days: t₅ = 5.39; P = 0.0029). In contrast, otolith-derived mm/d growth were between 13% and 79% lower than the laboratory growth model for both the last 8 days of growth and the penultimate 8 days of growth (paired t-test; last 8 days: t₅ = 6.54; P = 0.0012; penultimate 8 days: t₅ = 5.71; P = 0.0023). Across all years, field estimated growth rates were similar between heatwave conditions in 2014–2016 compared to more normal ocean conditions in 2010–2013 (1-way ANOVA; F₁,₁₀₂ = 0.75; P = 0.39), highlighting the difference in potential SGR predicted by the prSOM and the realized SGR observed in the field during heatwave conditions.

Figure 5. Comparison of juvenile Pacific Cod SGR (left panel;% wet weight/day) and mm/day growth rate (right panel) from individuals collected from Kodiak Island, AK in July 2010 and 2012–2016. Final 8 days of otolith growth prior to capture are represented by dark green triangles, and penultimate 8 days of otolith growth are represented by light green triangles. Laboratory models are represented by the blue growth curve and blue circles. Specific growth rate from the dominant stress-scapes in July near Kodiak Island are represented by gold diamonds. Error bars represent standard error. Gray polygons represent confidence and prediction intervals of the laboratory growth model.

Discussion

The application of seascape ecology methods to the evaluation of potential environmental stressors across a wide spatio-temporal scale on juvenile GOA Pacific Cod represents a new approach to visualize and synthesize the effects of extreme environmental variability on marine fisheries. Five dynamic stress-scapes were resolved from 2010 to 2016 that correspond to unique environmental conditions and growth responses of juvenile Pacific Cod. Stress-scape identity and spatial extent changed as a function of time, capturing unique conditions and responses during the 2014–2016 marine heatwave compared to the previous cooler conditions. Stress-scapes captured non-linear biological responses to the environment: when SGR was included as a variable (in addition to SST and pCO₂) there was a single stress-scape that dominated during heatwave years that did not dominate prior to heatwave years. This dominant stress-scape corresponds to observed elevated growth rates from the otolith data from Kodiak Island. Even though SGR is a function of SST, including fish growth as a variable in the stress-scape captures patterns of growth that are comparable to those derived from otoliths in situ.

Output from the stress-scape visualization indicated that juvenile Pacific Cod experienced moderately higher potential growth during heatwave years compared to more normal environmental conditions. Accelerated growth is a common bioenergetic response of teleost fishes to warmer ocean conditions, with increased temperatures contributing to increased metabolism and faster growth as long as caloric and oxygen requirements are adequately met, until an upper limit above which growth will begin to decline (e.g., Beauchamp, 2009; Thalmann et al., 2020). These results suggest that warmer temperatures during the juvenile life stage led to faster growth than would be expected during normal ocean conditions. However, in contrast with these observed potential growth rates, Pacific Cod abundance was anomalously low in the GOA during marine heatwave conditions. Pacific Cod, like many stenothermic fish species, are likely to experience life stage-specific thermal bottlenecks, especially in the embryonic and adult spawner stages where the fish may be particularly vulnerable to environmental variability (Dahlke et al., 2020). It has been shown that Pacific Cod hatch success is strongly correlated with temperature, with hatch success most optimal at temperatures between 3 and 6°C (Laurel and Rogers, 2020). During the 2014–2016 marine heatwave, optimal spawning habitat was restricted to shallower regions and earlier times of the year, and this decline in suitable spawning habitat may have led to widespread reductions in reproductive output. Those individuals that survived the thermal bottleneck at the embryonic life stage may have experienced higher growth rates during the larval and juvenile life stages in warmer conditions, especially if prey was available and “matched” in time and space (Laurel et al., 2021). Additionally, juvenile Pacific Cod collected in Kodiak nurseries were two to four months old by July and the product of selective mortality. In poor survival years, such as 2015, slower growing or smaller fish can experience higher rates of mortality, leaving only the faster growing fish in the population (Litvak and Leggett, 1992; Sogard, 1997). Thus, the results from the stress-scape visualization likely do not capture the full impact of the 2014–2016 marine heatwave on Pacific Cod simply because it captures only a snapshot of the entire Pacific Cod life history, but the approach highlights the spatial and temporal extent of anomalous and potentially stressful conditions. With the addition of information on prey quantity and quality and pCO₂ and temperature data throughout the water column, a more comprehensive exploration of the effects of the heatwave throughout the life history of Pacific Cod, and potentially other species, would be possible.

Otolith-derived SGR rates for field-collected juvenile Pacific Cod were between 39 and 187% higher than those predicted by the dominant stress-scape class for coastal Kodiak Island across most years of the analysis. This pattern highlights a potential underestimation of the stress-scape machine learning algorithm in predicting juvenile Pacific Cod growth in the field. The most parsimonious explanation for this result is that the temperature data used in the model did not accurately reflect shallow nursery conditions that often exceed 15°C in summer around Kodiak Island (Laurel et al., 2012). However, the laboratory model for SGR underestimated otolith-derived SGR across most years, while the laboratory model for mm/d growth overestimated otolith-derived growth. Growth in mass and growth in length are related but can reflect differences in energy allocation. For example, there is some evidence that seasonality can influence fish growth in the field, with fish adding more length earlier in the season when conditions are favorable for feeding, and then adding more mass later in the season to better prepare for winter conditions (e.g., Hurst and Conover, 2003; Hurst, 2007). Field growth may also exceed laboratory-based growth models due to a number of in situ processes, including higher food quality in the field (e.g., Mazur et al., 2007), favorable conditions in nearshore nursery habitat (e.g., Beck et al., 2001; Hinckley et al., 2019), and the effects of seasonality on growth rates, with juveniles feeding more vigorously in summer months to prepare for winter (e.g., Schwalme and Chouinard, 1999). In 2014, 2015, and 2012, two heatwave influenced years and a cool ocean year, respectively, juvenile Pacific Cod exhibited SGR in the field that were modestly higher than predicted from the laboratory temperature-dependent growth rate model. We note that the 2012 cohort of Pacific Cod was the strongest on record in the GOA since 1977, and growth was likely elevated in response to particularly favorable ocean conditions in that year and the production of multiple cohorts of juveniles (Hinckley et al., 2019). In the unfavorable ocean conditions occurring during the marine heatwave event in 2014 and 2015, only the fastest growing individuals were likely to survive to the juvenile life history stage, inflating the observed growth response of the survivors (e.g., Ricker, 1969; Litvak and Leggett, 1992; Moss et al., 2005). The variability between observed growth rates in the field and classified growth rates in the stress-scape suggests that future iterations of the stress-scape may benefit from added parameters related to selective mortality, foraging conditions, and additional environmental parameters related to ocean conditions.

In this analysis, we assume that the effects of temperature and pCO₂ on Pacific Cod growth are additive. We do not account for potential interaction between these two environmental variables as further studies are needed to elucidate the interactive effects between warming and ocean acidification on the growth of Pacific Cod early life stages. We acknowledge that the potential for interactive effects between temperature and pCO₂ on juvenile Pacific Cod growth may increase uncertainty in stress-scape classification and contribute to the deviation between Pacific Cod growth observed in the field compared to growth classified by the stress-scape. In general, early life stages of marine organisms tend to exhibit an increased sensitivity to the effects of ocean acidification in the presence of elevated water temperatures due to a heightened stress response (Harvey et al., 2013; Kroeker et al., 2013, 2014). However, fishes tend to be less sensitive than marine invertebrates to increased carbon dioxide in the water column even in the face of elevated temperature conditions due to their increased metabolism (e.g., Enzor et al., 2017), highly developed acid/base regulatory systems (e.g., Melzner et al., 2009), and potential for energy reallocation (e.g., Laubenstein et al., 2018). Studies examining the effects of elevated carbon dioxide concentrations on the early life stages of Alaskan gadids (including Walleye Pollock, a congener of Pacific Cod that shares a similar growth response at the juvenile stage), show limited impacts to growth response (Hurst et al., 2012, 2019). Therefore, we expect interactive effects of temperature and ocean acidification on SGR are likely minimal for juvenile Pacific Cod.

While previous dynamic classifications of the ocean include lower trophic level responses to environmental forcing (e.g., chl-a or carbon; Kavanaugh et al., 2014, 2018), the incorporation of fish physiology represents a notable revision and test of non-linear responses of higher trophic levels. Furthermore, our sensitivity analyses suggests that exclusion of biological responses results in an overly simple classification model. The bivariate prSOM (without growth rate) exhibited fewer output stress-scape classes compared to the full model and grouped several output classes together. In addition, the bivariate prSOM was unable to distinguish between heatwave conditions in 2014–2016 and more normal conditions in 2010–2013, with the dominant class remaining the same across all six years. Due to the quadratic relationship of Pacific Cod growth to ocean temperature (Hurst et al., 2010; Laurel et al., 2016), output classes in the full dynamic stress-scape model were able to successfully distinguish between similar growth rates that occurred in response to two different temperature ranges by classifying them into two separate output classes (classes 1 & 2). By including the biologically relevant growth response into the stress-scape, we could distinguish temperatures that were high enough to negatively impact growth rate from temperatures that were conducive to ideal growth. The utility of dynamic stress-scapes in capturing non-linear growth patterns could be expanded to include more complex biological patterns as indicators of changing ocean conditions, such as those associated with prey quality and quantity (e.g., Daly et al., 2017), predation (e.g., Fennie et al., 2020), or metabolic rate (e.g., Chung et al., 2019).

Machine learning models show promise in classifying marine heatwave events because they can generalize to regions of the ocean where there is limited data by being trained on more data-rich regions with similar conditions and then using the model to detect heatwave conditions. While there are several recent, well-received approaches to classifying marine heatwave events (e.g., Hobday et al., 2016, 2018; Oliver et al., 2018), many of these approaches exclude biological response in their classification schema for marine heatwaves. Our dynamic stress-scape classification of the GOA, which uses sea-surface measurements, successfully distinguished between marine heatwave conditions in 2014–2016 and more historically normal ocean conditions in 2010–2013. In general, stress-scape classes that corresponded to heatwave conditions exhibited higher SST, moderate pCO₂ concentrations, and elevated juvenile Pacific Cod SGR. Stress-scapes that corresponded to more normal ocean conditions exhibited moderate SST, lower pCO₂ concentrations, and moderate SGR. Dynamic stress-scapes effectively present the biological footprint of remotely sensed oceanographic data by characterizing spatio-temporal changes in the modeled SGR of juvenile Pacific Cod. The stress-scape framework shows promise in for merging spatio-temporal data and biological models to understand how heatwave conditions across a wide spatial and temporal range might affect living marine resources.

Our methods are not without limitations. The prSOM algorithm was chosen for this analysis because its topology preserving nature and its sensitivity to non-linear data. However, the prSOM implementation does not use k-fold cross validation to avoid overfitting, which would improve the classification results. Future efforts would include t-distributed stochastic neighbor embedding (t-SNE; Maatan and Hinton, 2008), which provides a well-defined “optimal solution,” lends itself to visualization of higher dimensional data (e.g., >3) and thus inclusion of additional in juvenile Pacific Cod fitness variables. Additionally, HAC can have issues when the geometry of the data is not flat (as is the case with the relationship between sea-surface temperature and growth rate). The DBSCAN (Ester et al., 1996) and OPTICS (Ankerst et al., 1999) algorithms could replace Hierarchical Agglomerative Clustering and be used in conjunction with a prSOM or t-SNE machine learning algorithm to provide a better stress-scape classification that distinguishes regions where the temperature is warm enough to constrain growth from regions where the temperature is conducive to optimal growth. Our visual results are a preliminary effort to introduce a visual system. In future work we will consider incorporating immersive visualization techniques for visual exploration of the ocean and coordinated multi-view visualization on the different types and sources of data for decision making. Efforts to expand on visualization within the fisheries community suggest interest in developing fast, scalable and explorative tools for fisheries scientists to understand the relationship between ocean variables, marine resources, and the algorithms used to process data from these two domains (Hermann and Moore, 2009).

Finally, our temperature-dependent growth model does not incorporate factors like prey quality, prey quantity, ocean currents, or ocean conditions at depth. Individual Based Modeling (IBM) and Regional Ocean Modeling Systems (ROMS) incorporate many of these factors and can be used to generate a spatial footprint of fish dispersal (Stockhausen et al., 2019). Information on the dominant dispersal pathways of Pacific Cod early life stages from an IBM/ROMS model could be paired with the stress-scapes by developing the duration of residence in each stress-scape. On a related note, IBM/ROMS models could be used to improve the data sampling used to train the prSOM. Currently, remotely sensed data are sampled uniformly at random at each epoch, where a different random sample is used at each iteration, but the distribution is uniform. The IBM/ROMS models could be used to estimate the probability of a juvenile Pacific Cod encountering potentially stressful conditions at each ocean pixel. This distribution could be used to inform the sampling method at each epoch of the prSOM so that the training data is more likely to represent the environmental context of the juvenile Pacific Cod. Our stress-scape framework is not intended as a substitute for an IBM, and it does not offer a comprehensive reconstruction of juvenile Pacific Cod growth history in the GOA (but see Hinckley et al. (2019) for an IBM currently in place for GOA Pacific Cod). Rather, the stress-scape represents a new way to visualize changes in both environmental conditions and Pacific Cod growing conditions before and during a marine heatwave event. Pacific Cod otoliths are used in the context of this study to support and serve as a point of comparison to the stress-scape classification of SST, pCO₂, and laboratory-derived juvenile Pacific Cod SGR.

In summary, stress-scapes offer three new and novel features relevant to the fields of marine resource management, ocean ecology and fisheries science. First, stress-scapes are based on both remotely sensed data and models that approximate biological response to remotely sensed environmental data. Second, machine learning models can generalize to regions of ocean where there is limited data by being trained on data that share similar conditions. Third, stress-scapes use machine learning algorithms that train on spatially relevant data, in this case based on temperature-growth models to approximate the growth rate of juvenile Pacific Cod. This produces an output that is potentially relevant for management in multiple fields and at multiple spatial and temporal scales. Environmentally, the stress-scapes are useful metrics that can benefit ecosystem-based and dynamic approaches to management of marine resources (Kavanaugh et al., 2016). From a management perspective, the biological relevance of the stress-scapes can be used to trace or even forecast potentially stressful conditions for a focus species. Downstream effects on socioeconomic conditions of resource users could also be identified by aligning the occurrence and longevity of stress-scape classes in space and time with economic indicators and extraction data. The process is extremely flexible and can handle large amounts of many types of time series data, furthering its utility.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

Ethical review and approval was not required for the animal study because the field work was conducted solely by NOAA’s Alaska Fisheries Science Center Laboratory in Newport, Oregon. Animal tissue, in this case otoliths, was provided to OSU several years after collection. This research was carried out in accordance with all applicable institutional and national guidelines at the time that the study was conducted; all work followed American Fisheries Society policies on the Guidelines for Use of Fishes in Research (https://fisheries.org/docs/policy_useoffishes.pdf) and AVMA (American Veterinary Medical Association) Guidelines on Euthanasia (olaw.nih.gov/sites/default/files/Euthanasia2007.pdf). There was no formal ethics review of this study because NOAA National Marine Fisheries Service does not have an Institutional Animal Care and Use Committee (IACUC) or an ethics approval processes for research on fishes and, at the time of these collection (2010–2016), OSU was not involved in the research.

Author Contributions

All authors conceived of and designed the project. JB, HT, and WK executed the research, analyzed the data, and wrote the manuscript with feedback from YZ, JM, BL, and MK. JB led the prSOM implementation, execution, analysis, and interpretation with input from YZ and MK. HT and BL compiled Pacific Cod growth data for incorporation into the prSOM. WK compiled environmental data for incorporation into the prSOM with assistance from MK. HT and JM completed otolith analyses and interpretation. BL provided fish samples.

Funding

This project was supported by the National Science Foundation Research Traineeship (NSF-NRT) award #1545188, Risk and Uncertainty Quantification and Communication in Marine Science and Policy. Research was also partially supported by the North Pacific Research Board grant no. 1903. JB, WK, and MK were supported by NASA awards 80NSSC18K0412 (Dynamic seascapes to support a biogeographic framework for a global marine biodiversity observing network) and 80NSSC20M0008 (Marine Biodiversity Observing Network in the Northern California Current: Understanding Patterns and Drivers of Biodiversity and Ecosystem Functioning from Plankton to Seascapes).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We greatly appreciated feedback on the development and implementation on this project from Lorenzo Ciannelli and the rest of the NSF-NRT leadership team at Oregon State University: Flaxen Conway, Ana Spalding, Julia Jones, and Katherine Hoffman. We are grateful to Tom Hurst for providing data and guidance related to Pacific Cod laboratory growth curves and to Thomas Murphy for assistance with processing the juvenile Pacific Cod otolith data used in this project. This effort was supported by a NSF-funded Research Traineeship at Oregon State University, NSF Grant 1545188 (NRT-DESE: Risk and Uncertainty Quantification in Marine Science). HT, JM, and BL were partially supported by the North Pacific Research Board grant no. 1903. In addition, JB, WK, and MK were supported by NASA awards 80NSSC18K0412 (Dynamic seascapes to support a biogeographic framework for a global marine biodiversity observing network) and 80NSSC20M0008 (Marine Biodiversity Observing Network in the Northern California Current: Understanding Patterns and Drivers of Biodiversity and Ecosystem Functioning from Plankton to Seascapes).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.656088/full#supplementary-material

Supplementary Video 1 is located here: https://figshare.com/articles/media/14356751 and Supplementary Video 2 is located here: https://figshare.com/articles/media/14356793.

References