In sean-rohan-NOAA/akfishcondition: Groundfish morphometric condition indicator

# Load packages
library(knitr)
library(akfishcondition)

region <- "BS"

make_map <- TRUE
pkg_version <- packageVersion("akfishcondition")

# Unzip map files
unzip(system.file("extdata/2022_ESR.zip", package = "akfishcondition"))

if(!dir.exists(here::here("plots", region))) {
  dir.create(here::here("plots", region), recursive = TRUE)
}

if(make_map){
  library(akgfmaps)

  ebs_layers <- akgfmaps::get_base_layers(select.region = "ebs",
                                          set.crs = "EPSG:3338")


  # Simple version
  survey_strata <- ebs_layers$survey.strata %>%
    dplyr::mutate(agg_stratum = Stratum) %>%
    dplyr::mutate(agg_stratum = replace(agg_stratum, agg_stratum %in% c(31,32), 30),
                  agg_stratum = replace(agg_stratum, agg_stratum %in% c(41,42,43), 40),
                  agg_stratum = replace(agg_stratum, agg_stratum %in% c(61,62), 60),
                  SURVEY = replace(SURVEY, SURVEY == "EBS_SHELF",  "EBS Shelf"), 
                  SURVEY = replace(SURVEY, SURVEY == "NBS_SHELF", "NBS")) %>% 
    dplyr::group_by(agg_stratum, SURVEY) %>% 
    dplyr::summarise()

  plot_ebs_nbs_survey_stations <- ggplot() +
    geom_sf(data = survey_strata,
            aes(fill = SURVEY),
            color = "black") +
    geom_sf_text(data = sf::st_centroid(survey_strata), 
                 aes(label = agg_stratum),
                 color = "black",
                 size = rel(6)) +
    geom_sf(data = sf::st_centroid(ebs_layers$survey.grid %>%
                                     dplyr::mutate(lab = "Station Location")),
            aes(shape = lab)) +
    ggplot2::geom_sf(data = ebs_layers$akland, 
                     fill = "grey85", 
                     color = "black") +
    geom_sf(data = ebs_layers$graticule,
            size = 0.2,
            alpha = 0.5) +
    ggplot2::coord_sf(xlim = ebs_layers$plot.boundary$x, 
                      ylim = ebs_layers$plot.boundary$y) +
    ggplot2::scale_x_continuous(name = "Longitude", 
                                breaks = ebs_layers$lon.breaks) + 
    ggplot2::scale_y_continuous(name = "Latitude", 
                                breaks = ebs_layers$lat.breaks) +
    scale_shape_manual(values = "x") +
    scale_fill_manual(values = c("#0085CA", "#54ADDB")) +
    theme_bw() +
    ggplot2::theme(axis.title = element_blank(),
                   axis.text = element_text(color = "black"),
                   axis.ticks = element_line(color = "black"),
                   panel.border = element_rect(color = "black", fill = NA),
                   panel.grid = element_blank(),
                   panel.background = element_rect(color = "black", fill = "white"),
                   legend.title = element_blank(),
                   legend.margin = margin(-12,0,0,0),
                   legend.position = c(0.2,0.1),
                   legend.background = element_blank())

  ragg::agg_png(filename = here::here("plots", region, "EBS_NBS_survey_area.png"), width = 5, height = 5, units = "in", res = 600)
  print(plot_ebs_nbs_survey_stations)
  dev.off()
}

Contributed by Bianca Prohaska^1^, Rebecca Howard^1^ and Sean Rohan^1^

^1^Resource Assessment and Conservation Engineering Division, Alaska Fisheries Science Center, National Marine Fisheries Service, NOAA
Contact: sean.rohan@noaa.gov Last updated: September 2024

Description of Indicator: Length-weight residuals represent how heavy a fish is per unit body length and are an indicator of somatic growth variability [@Brodeur2004]. Therefore, length-weight residuals can be considered indicators of prey availability, growth, general health, and habitat condition [@Blackwell2000; @Froese2006]. Positive length-weight residuals indicate better condition (i.e., heavier per unit length) and negative residuals indicate poorer condition (i.e., lighter per unit length) [@Froese2006]. Fish condition calculated in this way reflects realized outcomes of intrinsic and extrinsic processes that affect fish growth, which can have implications for biological productivity through direct effects on growth and indirect effects on demographic processes, such as reproduction and mortality (e.g., [@Rodgveller2019; @Barbeaux2020]).

```rFigure 1. AFSC/RACE GAP summer bottom trawl survey strata (10-90) and station locations (x) in the eastern Bering Sea and northern Bering Sea. ", echo=FALSE} knitr::include_graphics(path = here::here("plots", region, "ebs_nbs_survey_area.png"))

The groundfish morphometric condition indicator is calculated from paired fork lengths (mm) and weights (g) of individual fish that were collected during bottom trawl surveys of the eastern Bering Sea (EBS) shelf and northern Bering Sea (NBS), which were conducted by the Alaska Fisheries Science Center’s Resource Assessment and Conservation Engineering (AFSC/RACE) Groundfish Assessment Program (GAP). Fish condition analyses were applied to walleye pollock (_Gadus chalcogrammus_), Pacific cod (_Gadus macrocephalus_), arrowtooth flounder (_Atheresthes stomias_), yellowfin sole (_Limanda aspera_), flathead sole (_Hippoglossoides elassodon_), northern rock sole (_Lepidopsetta polyxystra_), and Alaska plaice (_Pleuronectes quadrituberculatus_) collected in bottom trawls at standard survey stations (Figure 1). Alaska plaice were not sampled in 2024. For these analyses and results, survey strata 31 and 32 were combined as stratum 30; strata 41, 42, and 43 were combined as stratum 40; and strata 61 and 62 were combined as stratum 60. Northwest survey strata 82 and 90 were excluded from these analyses.

To calculate indicators, length-weight relationships were estimated from linear regression models based on a log-transformation of the exponential growth relationship, _W_ = _aL_^_b_^, where _W_ is weight (g) and _L_ is fork length (mm) for all areas for the period 1997–2024 (EBS: 1997–2024, NBS: 2010, 2017, 2019, 2021-2023). Unique  intercepts (*a*) and slopes (*b*) were estimated for each survey stratum, sex, and interaction between stratum and sex to account for sexual dimorphism and spatial-temporal variation in growth and bottom trawl survey sampling. Length-weight relationships for 100–250 mm fork length walleye pollock (corresponding with ages 1–2 years) were calculated separately from adult walleye pollock (> 250 mm). Residuals for individual fish were obtained by subtracting observed weights from bias-corrected weights-at-length that were estimated from regression models. Length-weight residuals from each stratum were aggregated and weighted proportionally to total biomass in each stratum from area-swept expansion of mean bottom-trawl survey catch per unit effort (CPUE; i.e., design-based stratum biomass estimates). Variation in fish condition was evaluated by comparing average length-weight residuals among years. To minimize the influence of unrepresentative samples on indicator calculations, combinations of species, stratum, and year with a sample size <10 were used to fit length-weight regressions, but were excluded from calculating length-weight residuals for both the EBS and NBS. Morphometric condition indicator time series, code for calculating the indicators, and figures showing results for individual species are available through the *akfishcondition* R package and GitHub repository (https://github.com/afsc-gap-products/akfishcondition).

**Status and Trends**: Fish condition, indicated by length-weight residuals, has varied over time for all species examined in the EBS (Figure 2 & 3). Since 2022, a declining trend in condition was observed in small walleye pollock (100-250 mm), northern rock sole, and Pacific cod; however, these species were still within one standard deviation of the historical mean (Figure 2). Large walleye pollock (>250 mm), yellowfin sole, and arrowtooth flounder exhibited an increase in condition from 2023 to 2024, while little to no change in condition was observed in flathead sole (Figure 2). Despite observing increasing or decreasing trends in condition among most of the species examined, the 2024 mean condition of all species, except for northern rock sole, fell within two standard errors of the mean condition in 2023 (Figure 2). Furthermore, the mean condition of every species fell within one standard deviation of the historical mean in 2024 (Figure 2).

In the EBS in 2024, condition was negative for all species examined across most strata, except for yellowfin sole on the inner and middle shelf (Strata 10 and 30), and arrowtooth flounder on the southern outer shelf (Stratum 50) (Figure 3). 

```rFigure 2. Morphometric condition of groundfish species collected during AFSC/RACE GAP standard summer bottom trawl surveys of the eastern Bering Sea shelf (1999 to 2024) based on residuals of length-weight regressions. The dash in the blue boxes denotes the mean for that year, the box denotes one standard error, and the lines on the boxes denote two standard errors. Lines on each plot represent the historical mean, dashed lines denote one standard deviation, and dotted lines denote two standard deviations.",  message=FALSE, warning=FALSE}
# Set factor levels for plotting order
fig2  <- plot_anomaly_timeseries(x = EBS_INDICATOR$FULL_REGION,
                        region = "BS",
                        fill_color = "#0085CA",
                        var_y_name = "mean_wt_resid",
                        var_y_se_name = "se_wt_resid",
                        var_x_name = "year",
                        y_title = "Length-weight residual (ln(g))",
                        plot_type = "box",
                        format_for = "rmd")

fig2_png  <- plot_anomaly_timeseries(x = EBS_INDICATOR$FULL_REGION,
                        region = "BS",
                        fill_color = "#0085CA",
                        var_y_name = "mean_wt_resid",
                        var_y_se_name = "se_wt_resid",
                        var_x_name = "year",
                        y_title = "Length-weight residual (ln(g))",
                        plot_type = "box",
                        format_for = "png")

print(fig2 + theme_condition_index())

png(here::here("plots", region, "EBS_condition.png"),width=6,height=7,units="in",res=600)
print(fig2_png + theme_blue_strip())
dev.off()

```rFigure 3. Length-weight residuals by survey stratum (10-60) for seven eastern Bering Sea shelf groundfish species and age 1–2 walleye pollock (100–250 mm) sampled in the AFSC/RACE GAP standard summer bottom trawl survey, 1999-2024. Length-weight residuals are not weighted by stratum biomass.", message=FALSE, warning=FALSE}

EBS stratum-biomass weighted residual stratum stacked bar plots

fig3 <- akfishcondition::plot_stratum_stacked_bar(x = EBS_INDICATOR$STRATUM, region = region, var_x_name = "year", var_y_name = "stratum_resid_mean", y_title = "Length-weight residual (ln(g))", var_group_name = "stratum", fill_title = "Stratum", fill_palette = "BrBG")

fig3_png <- akfishcondition::plot_stratum_stacked_bar(x = EBS_INDICATOR$STRATUM, region = region, var_x_name = "year", var_y_name = "stratum_resid_mean", y_title = "Length-weight residual (ln(g))", var_group_name = "stratum", fill_title = "Stratum", fill_palette = "BrBG")

print(fig3 + theme_condition_index()+ theme(strip.text = element_text(size = 20)))

```r
png(here::here("plots", region, "EBS_SBW_stratum_stacked_bar.png"), width = 6, height = 7, units = "in", res = 600)
print(fig3_png + theme_blue_strip() + theme(legend.position = "right", legend.title = element_text(size = 9)))
dev.off()

The NBS was not surveyed in 2024; however, in the NBS in 2023, positive condition was observed for large walleye pollock (>250 mm), which has been increasing since 2021. The remaining species exhibited near-average condition in the NBS in 2023, except for yellowfin sole, which exhibited negative condition and has been declining since 2019 (Figure 4).

In 2023 large walleye pollock (>250 mm) condition was positive in all NBS strata, whereas condition was previously negative in all strata from 2021-2022 (Figure 5). Pacific cod, small walleye pollock (100-250 mm), Alaska plaice, and yellowfin sole condition have been consistently negative across all strata since 2021, with a notable exception in 2023 of positive condition for Pacific cod in the inner southern NBS shelf, and Alaska plaice in the northern inner NBS shelf and Norton Sound (Figure 5).

```rFigure 4. Morphometric condition of groundfish species collected during AFSC/RACE GAP standard summer bottom trawl surveys of the northern Bering Sea shelf (2010, 2017, 2019 and 2021-2023) based on residuals of length-weight regressions. The dash in the blue boxes denotes the mean for that year, the box denotes one standard error, and the lines on the boxes denote two standard errors. Lines on each plot represent the historical mean, dashed lines denote one standard deviation, and dotted lines denote two standard deviations.", message=FALSE, warning=FALSE}

Set factor levels for plotting order

fig4 <- plot_anomaly_timeseries(x = NBS_INDICATOR$FULL_REGION, region = "BS", fill_color = "#54ADDB", var_y_name = "mean_wt_resid", var_y_se_name = "se_wt_resid", var_x_name = "year", y_title = "Length-weight residual (ln(g))", plot_type = "box", format_for = "rmd")

fig4_png <- plot_anomaly_timeseries(x = NBS_INDICATOR$FULL_REGION, region = "BS", fill_color = "#54ADDB", var_y_name = "mean_wt_resid", var_y_se_name = "se_wt_resid", var_x_name = "year", y_title = "Length-weight residual (ln(g))", plot_type = "box", format_for = "png")

print(fig4 + theme_condition_index())

```r
png(here::here("plots", region, "NBS_condition.png"), width=6, height=7, units="in", res=600)
print(fig4_png + theme_blue_strip())
dev.off()

```rFigure 5. Length-weight residuals by survey stratum (70, 71 and 81) for four northern Bering Sea shelf groundfish species and age 1–2 walleye pollock (100–250 mm) sampled in the AFSC/RACE GAP standard summer bottom trawl survey during 2010, 2017, 2019 and 2021-2023. Length-weight residuals are not weighted by stratum biomass.", message=FALSE, warning=FALSE}

NBS stratum-biomass weighted residual stratum stacked bar plots

fig5 <- akfishcondition::plot_stratum_stacked_bar(x = NBS_INDICATOR$STRATUM, region = region, var_x_name = "year", var_y_name = "stratum_resid_mean", y_title = "Length-weight residual (ln(g))", var_group_name = "stratum", fill_title = "Stratum", fill_palette = "PuBu")

fig5_png <- akfishcondition::plot_stratum_stacked_bar(x = NBS_INDICATOR$STRATUM, region = region, var_x_name = "year", var_y_name = "stratum_resid_mean", y_title = "Length-weight residual (ln(g))", var_group_name = "stratum", fill_title = "Stratum", fill_palette = "PuBu")

print(fig5 + theme_condition_index() + theme(strip.text = element_text(size = 20)))

```r
png(here::here("plots", region, "NBS_SBW_stratum_stacked_bar.png"), width = 6, height = 7, units = "in", res = 600)
print(fig5_png + theme_blue_strip() + theme(legend.position = "right", legend.title = element_text(size = 9)))
dev.off()

Factors influencing observed trends: Similar to trends observed in 2023, temperature appears to influence morphological condition of species in the EBS, so near-average cold pool extent and water temperatures in 2024 likely played a role in the near-average condition (within 1 S.D. of the mean) for all species examined in the EBS. Historically, particularly cold years tend to correspond with negative condition, while particularly warm years tend to correspond to positive condition. For example, water temperatures were particularly cold during the 1999 Bering Sea survey, a year in which negative condition was observed for all species for which data were available. In addition, spatiotemporal factor analyses suggest the morphometric condition of age-7 walleye pollock is strongly correlated with cold pool extent in the EBS [@Gruss2021]. In recent years, warm temperatures across the Bering Sea shelf, since the record low seasonal sea ice extent in 2017–2018 and historical cold pool area minimum in 2018 [@Stabeno2019a], may have influenced the positive trend in the condition of several species from 2016 to 2019. However, despite near-average temperature in 2023 large walleye pollock (>250 mm) condition in the EBS was the second lowest recorded over the time series. In 2024 temperature was near-average again, and while large walleye pollock (>250 mm) condition was still negative, it was less negative than that observed in 2023.

Although warmer temperatures may increase growth rates if there is adequate prey to offset temperature-dependent increases in metabolic demand, growth rates may also decline if prey resources are insufficient to offset temperature-dependent increases in metabolic demand. The influence of temperature on growth rates depends on the physiology of predator species, prey availability, and the adaptive capacity of predators to respond to environmental change through migration, changes in behavior, and acclimatization. For example, elevated temperatures during the 2014–2016 marine heatwave in the Gulf of Alaska led to lower growth rates of Pacific cod and lower condition because available prey resources did not make up for increased metabolic demand [@Barbeaux2020].

Other factors that could affect morphological condition include survey timing, stomach fullness, fish movement patterns, sex, and environmental conditions [@Froese2006]. The starting date of annual length-weight data collections has varied from late May to early June and ended in late July-early August in the EBS, and mid-August in the NBS. Although we account for some of this variation by using spatially-varying coefficients in the length-weight relationship, variation in condition could relate to variation in the timing of sample collection within survey strata. Survey timing can be further compounded by seasonal fluctuations in reproductive condition with the buildup and depletion of energy stores [@Wuenschel2019]. Another consideration is that fish weights sampled at sea include gut content weights, so variation in gut fullness may influence weight measurements. Since feeding conditions vary over space and time, prey consumption rates and the proportion of total body weight attributable to gut contents may be an important factor influencing the length-weight residuals.

Finally, although the condition indicators characterize temporal variation in morphometric condition for important fish species in the EBS and NBS, they do not inform the mechanisms or processes behind the observed patterns.

Implications: Fish morphometric condition can be considered an indicator of ecosystem productivity with implications for fish survival, maturity, and reproduction. For example, in Prince William Sound, the pre-winter condition of herring may determine their overwinter survival [@Paul1999], differences in feeding conditions have been linked to differences in morphometric condition of pink salmon in Prince William Sound [@Boldt2004], variation in morphometric condition has been linked to variation in maturity of sablefish [@Rodgveller2019], and lower morphometric condition of Pacific cod was associated with higher mortality and lower growth rates during the 2014–2016 marine heat wave in the Gulf of Alaska [@Barbeaux2020]. Condition can also be an indicator of stock status relative to carrying capacity because morphometric condition is expected to be high when the stock is at low abundance and low when the stock is at high abundance because of the effects of density-dependent competition [@Haberle2023]. Thus, the condition of EBS and NBS groundfishes may provide insight into ecosystem productivity as well as fish survival, demographic status, and population health. However, survivorship is likely affected by many factors not examined here.

Another important consideration is that fish condition was computed for all sizes of fishes combined, except in the case of walleye pollock. Examining condition of early juvenile stage fishes not yet recruited to the fishery, or the condition of adult fishes separately, could provide greater insight into the value of length-weight residuals as an indicator of individual health or survivorship [@Froese2006], particularly since juvenile and adult walleye pollock exhibited opposite trends in condition in the EBS this year.

The near-average condition for most species in 2024 may be related to the near historical average temperatures observed. However, trends in recent years such as prolonged warmer water temperatures following the marine heat wave of 2014-16 [@Bond2015] and reduced sea ice and cold pool areal extent in the eastern Bering Sea [@Stabeno2019a] may affect fish condition in ways that have not yet been determined. Additionally, periods of high fishing mortality that reduce population biomass are likely to increase body condition because of the compensatory alleviation of density-dependent competition [@Haberle2023]. As we continue to add years of length-weight data and expand our knowledge of relationships between condition, growth, production, survival, and the ecosystem, these data may increase our understanding of the health of fish populations in the EBS and NBS.

Research priorities: Recent studies have sought to improve the understanding of factors that affect fish growth and morphometric condition in the EBS (e.g.[@Oke2022]) although significant knowledge gaps remain for numerous species and life history stages. Research is underway to evaluate connections between morphometric condition indices, temperature, and density-dependent competition.