Introduction

External pigmentation in fishes is under evolutionary pressure and has consequences for fitness, mediated through effects of coloration on survival and reproduction (Leclercq et al. 2010). While fish also vary in their internal (flesh) coloration, less is known about the mechanisms and selective pressures influencing pigmentation under the skin. Flesh coloration in fish is known to be affected by both intrinsic (hormones and nervous system) and extrinsic factors (chemicals and food; Prince 1916; Tunison et al. 1947). Color production in fishes is generally either due to pigmentation or is structurally based (Price et al. 2008), the latter being more common in blue-colored fishes (Fujii 1993). Variation in flesh color has been most extensively studied for species that are raised in aquaculture settings, where deviation from expected color will affect the market value; for example, salmon are known to vary in their flesh redness both as a function of diet and genetics (Torrissen and Naevdal 1988).

Blue coloration in tissues and serum has been noted in a range of fishes (reviewed in Nomura 1967), both in freshwater and marine systems. In sculpins (Cottidae), it has been shown that blue coloration is due to the bile pigment biliverdin, a byproduct of heme catabolism; the coloration can be found in the serum (Bada 1970; Low and Bada 1974) and the tissues (Yamaguchi et al. 1976). Bard et al. (1996) described an unusual series of “green” tunas (Scombridae; multiple genera) landed in ports of West Africa, which had confirmed biliverdin in blue-stained bone structure. Some eels (Anguilla japonica and A. rostrata; Yamaguchi et al. 1966; Ellis and Poluhowich 1981) and wrasses (Labrus bergylta; Clark et al. 2016) are known to have blue–green color serum and high biliverdin concentrations in plasma. Blue flesh (i.e., muscle tissue) has also been reported in the humphead wrasse (Cheilinus undulatus; Yamaguchi and Matsura 1969), garfish (Belone belone; Fontaine 1941), and Pacific saury (Cololabis saira; Nomura 1967). Blue-morphs of lake-dwelling walleye (Sander vitreus) are caused by biliverdin-bound fluorescent proteins in the mucus layer of the fish, which has been suggested as a response to UV radiation (Schaefer et al. 2014; Ghosh et al. 2016). Blue blood and tissues, as a result of high biliverdin concentrations, has been documented in other vertebrates, including lizards (Prasinohaema and other genera Rodriguez et al. 2018), frogs (Chiromantis; Grismer et al. 2007), and worms (Manduca sexta; Goodman et al. 1985).

Fishers and ecologists alike have long puzzled over the striking bright blue or blue–green coloration within the tissues of some benthic fishes in the greenling (Hexagrammidae) and sculpin families (Yamaguchi et al. 1976; Love 2011). In several of these species that support recreational and commercial fisheries in the northeastern Pacific—lingcod (Ophiodon elongatus), rock greenling (Hexagrammos lagocephalus), and cabezon (Scorpaenichthys marmoratus)—the coloration of the flesh is at times bright enough that it can be observed on the outside of the fish (Fig. 1). In addition, blue-fleshed fish can be particularly valued in the marketplace (Chuck’s Seafood, Charleston, Oregon, personal communication to AWEG; Bland 2014). However, there are no published quantitative analyses that examine incidence, and proximate or ultimate factors that cause blue coloration in the greenling family. Lingcod provide a particularly valuable opportunity for understanding the ecological factors explaining the prevalence of blue flesh in fish across a large geographic range, due to their cosmopolitan distribution throughout the northeastern Pacific Ocean, from the Gulf of Alaska to central Baja Mexico (Miller and Lea 1972).

Fig. 1
figure 1

Examples of blue-fleshed fish and serum. a Blue-fleshed rock greenling (Hexagrammos lagocephalus; Aleutian Islands, AK; Photo by Scott Gabara). b, c Blue flesh and mouth tissues of a lingcod (Port Orford, OR; Photos by Aaron Galloway); d Blue serum from lingcod (San Juan Islands, WA; Photo by Aaron Dufault); e Representative range of color variants—from the left, individuals 1, 3, and 4 were designated ‘brown’ and individuals 2 and 5 were designated ‘blue’ (Photo by Laurel Lam)

Lingcod reside in diverse habitats from the shoreline to a depth of at least 427 m (Miller and Geibel 1973), and typically exhibit a high degree of site fidelity (Starr et al. 2004; Greenley et al. 2016). Lingcod make a seasonal migration from deeper to shallower nearshore waters where they spawn in late fall through winter (Phillips 1958; Starr et al. 2004). As lingcod age and grow, they tend to move to deeper, moderate- to high-relief rocky habitats (Cass et al. 1990; Bassett et al. 2018). Females grow faster than males and reach larger sizes, usually inhabiting deeper depths (Miller and Geibel 1973). Lingcod are generalist predators that consume a diversity of pelagic and demersal fishes and mobile invertebrates (Beaudreau and Essington 2007, 2009), including cephalopods and crustaceans that contain respiratory proteins (hemocyanins) that are blue when oxygenated (van Holde et al. 2001). Blue flesh in lingcod is a frequently discussed topic in the fishing community and is often attributed to differences in the diets of individual fish; however, this assertion has not, to our knowledge, been tested formally.

A valuable method for understanding diets in marine fishes is through quantification of the concentrations and composition of fatty acids in a consumer’s tissues (Parrish 2013; Pethybridge et al. 2018). Individual fatty acids and key classes of physiologically important fatty acids are commonly used as trophic tracers (i.e., ‘biomarkers’) for inferring feeding relationships among predators and prey (Budge et al. 2002; Dalsgaard et al. 2003). The integration of fatty acids into the tissues of fishes is complex and depends largely upon their diet (Xu et al. 2020) (exogenous sources) but also depends on individual species physiology (endogenous processes), which is also governed by environmental conditions (reviewed in Arts and Kohler 2009; Parrish 2013; Xu et al. 2020). Essential fatty acids (EFAs) such as docosahexaenoic (DHA; 22:6ω-3), eicosapentaenoic (EPA; 20:5ω-3), and arachidonic (ARA; 20:4ω-6) fatty acids are important for the nutritional physiology of fishes (Sargent et al. 1999; Copeman et al. 2002; Bell and Sargent 2003; Parrish 2013). Moreover, the ratio of ω-3:ω-6 fatty acids, and total fatty acids can be used as indices of fish condition (Pethybridge et al. 2015). A common application of fatty acids as trophic indicators in fishes is the comparative analysis of conspecifics between life stages, with different coloration, or based on differences in the collection habitats, or across spatial or environmental gradients (Budge et al. 2002; Pethybridge et al. 2015; Copeman et al. 2016). This kind of comparative approach is often used for species in which consumer-specific fatty acid trophic modification is unknown (Galloway and Budge 2020), and does not aim to infer consumer diet composition, but rather assumes that differences in fatty acids of conspecifics between treatments or across a geographical gradient indicates dissimilarities in trophic ecology.

We used data from multiple field collection efforts to broadly document spatial and demographic patterns in the prevalence of blue coloration in lingcod throughout their range along the U.S. west coast and to investigate potential physiological and environmental correlates of blue flesh in the species. We performed an initial pilot study based on lingcod caught in southern Oregon in 2016 and tested whether the fatty acids differed between blue and non-blue (hereafter, ‘brown’) fish. We then extracted the fatty acids from a subset of the lingcod sampled in 2015–2017 over a broad geographic range, from southeastern Alaska to southern California, to determine whether blue and brown lingcod could be distinguished based on their fatty acid profiles. We used the full latitudinal dataset to evaluate the relative importance of influential fatty acids (as determined by the trophic tracer analyses), geographic factors, and biological variables in determining the probability of blue flesh in lingcod.

Methods

Field and laboratory methods

Southern Oregon survey

We collected muscle tissue, using an 8 mm biopsy punch, from the left dorsal area of lingcod (n = 72) in the Redfish Rocks Marine Reserve and Humbug Comparison Area near Port Orford, Oregon. The samples were collected by the Oregon Department of Fish and Wildlife Marine Reserves Program on six dates in June and July 2016 using longline and hook-and-line sampling methods described in Huntington and Watson (2017). In the field, fork length (cm), sex (male/female), muscle tissue coloration (blue, brown), fish condition comments, bottom depth, and collection site coordinates were recorded. We stored all fish tissues on ice until they could be frozen at the lab (e.g., < 6 h after collection) in an − 80 °C freezer. This dataset was used to identify emergent patterns in the sex ratio, fatty acid composition, and comparisons of total ω-6 fatty acids (see methods below) between blue and brown fish. This initial analysis thus identified hypotheses to test with the independent and much larger fatty acid dataset from the latitudinal survey.

Latitudinal survey

Lingcod were collected in 2015–2017 from chartered commercial passenger fishing vessels from Yakutat, Alaska to San Diego, California, USA (54°30’ N–34°30’ N). The following seven regions were broadly identified (Fig. 2): southeast Alaska (54°30’ N–59°48’ N), Puget Sound, coastal Washington (46°16’ N–49° N), Oregon (42° N–46°16’ N), northern California (38°02’ N–42° N), central California (34°30’ N–38°02’ N), and southern California (32°32’ N–34°30’ N). Three to four fishing ports were selected per region (24 ports total) based on location and accessibility of charter vessels. Regions were generally selected based on major biogeographic features (point conception, Cape Mendocino, Columbia River, and Puget Sound). We collected lingcod using hook-and-line (jigging with artificial lures). Hook-and-line gear is an effective method of capturing fishes in high-relief, rocky habitat that cannot be sampled using nets (e.g., trawls) and is commonly used in lingcod studies (e.g., Beaudreau and Essington 2007, 2009). Like all gear types, hook-and-line can potentially introduce selection bias, such as the potential to capture more actively feeding fish and selecting for post-settlement juvenile and adult lingcod (> 20 cm total length), as smaller lingcod have not yet settled into rocky reef habitats. Acknowledging these limitations, we aimed to capture a representative sample of lingcod across a wide range of size, sex, and age classes by sampling shallow (< 61 m) and deep (61–168 m) nearshore and offshore rocky reefs equally. We fished for an average of 2.5 days out of each port to obtain the targeted sample size of 75–100 lingcod per port. Additional samples were provided by the Alaska commercial longline fishery, a rockfish bycatch study in Puget Sound (K. Andrews, unpublished data), the Oregon Department of Fish and Wildlife Marine Reserves Program (Huntington et al. 2014), and the California Collaborative Fisheries Research Program (Wendt and Starr 2009).

Fig. 2
figure 2

Map summarizing the regional proportions of blue and brown lingcod by sex. The percentage of all fish that were blue is shown for each region in parentheses. Additional details about the sample sizes in each category and numbers included in the fatty acids analysis are reported in Table S1

Location information (i.e., latitude, longitude, bottom depth, relief, and time of day) was collected at the start and stop of each fishing drift. Total length (cm), weight (kg), sex, and color were recorded for each individual immediately upon capture. Samples were then placed on ice or frozen until dissection. Sex was determined externally (by the presence or absence of a conical papilla) and internally (by the presence or absence of ovaries or testes). Color was visually assessed both externally and internally by two authors to ensure consistency (BLB, LSL). We determined external color based on skin pigmentation immediately after each fish was landed, then verified external color by examining the buccal cavity and dorsal fin membrane color. This secondary check was done because mottling patterns on the skin can confound paler blue hues, while the buccal cavity and fin membrane are typically uniform in color and offer a cleaner presentation of the blue coloration (Fig. 1). Internal color provided a third check on coloration and was determined based on muscle tissue during the time of dissection, between 8 and 36 h after capture. In externally blue fish, the internal color of the flesh is also blue; however, in the more common externally brown fish, the internal color of the flesh is white/gray. While the intensity of blue pigmentation was observed as a gradient, color was recorded as binary (blue vs. brown; Fig. 1e) to avoid subjective determination of color intensity. The multiple methods (skin, buccal cavity, fin membrane, and flesh) used to assess and verify color improved confidence in our designation of external and internal color. During the full dissection process, we recorded liver, gonad, whole stomach weight, and stomach content weight (g). A piece of muscle tissue was collected from the anterior dorsal region on the left side for each individual and placed in a 1.8 mL cryogenic vial and frozen at  –20 °C until fatty acids extraction in the spring of 2018. A total of 2157 fish with color, sex, and size information were recorded.

Lab analysis of fatty acids

As it was not possible to extract fatty acids from all latitudinal survey lingcod tissue samples, we selected a subset of 175 samples for fatty acid analysis (Table S1) using a stratified random sampling approach. Since an initial analysis of the southern Oregon data indicated that blue coloration varied by sex, we subsampled to achieve approximately equal proportions of sex and coloration (blue or brown) within each of the seven regions sampled in the latitudinal survey. A minimum of 5 lingcod were sampled from each stratum (i.e., combination of geographical region, sex, and flesh color), except for the case of blue male lingcod in Washington, where we did not have any representative fish tissues to draw from (Table S1). All samples were lyophilized, while frozen, for 48 h, to achieve complete dryness. Total fatty acid extraction procedures followed Thomas et al. (2020), using a ~ 20 mg dry weight tissue aliquot. Briefly, all lipids were extracted from fish tissues using a 2:1:0.75 mixture of chloroform:methanol:NaCL (0.9%). An unmethylated internal standard (19:0) was added to each sample to quantify extraction and methylation efficiency. Samples were sonicated to facilitate cell lysis and centrifuged to isolate the lipid-containing organic phase. An aliquot of the organic phase was removed for gravimetric determination of total lipid content (mg lipid g−1 dry wt.) following Kainz et al. (2017). The remaining portion of the organic phase was evaporated to dryness and 1 mL toluene and 2 mL methanolic sulfuric acid (1%) were added immediately before incubating samples at 90 °C for 90 min to form fatty acid methyl esters. After incubation, purified samples were neutralized with 1.5 mL KHCO3 (2%) and extracted twice with n-hexane. Analysis of fatty acid methyl esters followed procedures described in Taipale et al. (2013) using a Shimadzu QP2020 gas chromatograph-mass spectrometer at the Oregon Institute of Marine Biology.

Data analysis

We used multiple datasets and analytic methods to address the study objectives. A visual summary of datasets used, study questions, statistical methods, and key results is shown in Figure S1.

Data preparation

Prior to analyses, we prepared the data as follows. Out of the 2157 lingcod from the latitudinal survey, there was agreement between external and internal color designations for 1934 fish (89.7%). There were discrepancies between external and internal color for 136 fish (6.3%); specifically, 2.4% of the fish were assigned brown externally and blue internally, whereas 3.9% of the fish were assigned blue externally and brown internally. These 136 fish were excluded from further analysis. Twenty-two fish were missing an external color designation (1.0%) and 65 fish were missing an internal color designation (3.0%), due to accidental omission while recording field data. Rather than reducing the sample by 87 fish, we assigned the 22 fish of unknown external color to their observed internal color. All further analyses were performed using external color designations. A total of 2021 lingcod with complete size, sex, and color information were available for analysis.

Fatty acids

A total of 46 unique fatty acids were identified and used in our analyses of both fatty acid datasets. For the Southern Oregon survey, we compared the sex ratio of blue and brown fish using a Chi-square test and analyzed fatty acids of 72 lingcod (41 male and 31 female) in 2017. These preliminary analyses inspired the larger and more robust latitudinal lingcod fatty acids effort, designed to account for blueness, sex, and region (described below). We compared multivariate fatty acid composition using two-way permutational multivariate analysis of variance (PERMANOVA; 9999 permutations, Type III SS, p values calculated from unrestricted permutation of the raw data; Anderson et al. 2008) on the factors ‘color’ and ‘sex’. Concentrations of key a-priori determined fatty acid categories [total polyunsaturated fatty acids (PUFA), total ω-6 fatty acids, the ω-3:ω-6 fatty acid ratio, total long-chain PUFA [LCPUFA; sum of 20:4ω-6 (ARA), 20:5ω-3 (EPA), 22:6ω-3 (DHA)], and total lipids] were log-transformed and compared with two-tailed t tests.

To further examine differences in fatty acid composition and concentrations at a larger geographic scale (i.e., combinations of region, sex, and color), we extracted and analyzed 175 latitudinal fatty acid samples. Two of these samples were later excluded due to their extreme outlier status when visualized in multivariate space; on further inspection, these samples had anomalous and unreasonable levels of several fatty acids, likely due to extraction process errors. We, therefore, performed subsequent multivariate analyses on 173 samples using 3-way PERMANOVA, using log-transformed fatty acid concentration data {log[X + 1]}, with the factors ‘region’ (random factor), ‘sex’, and ‘color’ (both fixed factors). We focused analyses on the fatty acid concentration data because concentrations are more informative from a physiological standpoint. Table S2 reports the mean and SD of all fatty acids summarized by lingcod flesh color. We tested the homogeneity of multivariate dispersions (PERMDISP) of factors analyzed in the PERMANOVA. We used multivariate similarity percentages (SIMPER) on the transformed fatty acid concentration data to identify the fatty acids that were most influential in separating the blue and brown fish. In tables, we report the untransformed, actual mean concentrations of fatty acids for ease of interpretation, expressed as mg fatty acids g−1 dry weight. We visualized the multivariate group dispersion patterns using a non-metric multidimensional scaling (nMDS) analysis (Euclidean distance), with individual fatty acids correlated with the ordination axes shown as vector overlays. Individual and summary categories of fatty acids that were identified in the initial southern Oregon regional collections, in the literature (e.g., Pethybridge et al. 2015), or the SIMPER analysis on the full west coast latitudinal collections were compared with two-sample t tests and boxplots. We used Primer v6 (Clarke and Gorley 2006) with the PERMANOVA + add-on (Anderson et al. 2008) for multivariate analyses. All tests of statistical significance were two-tailed and generally used α = 0.05. For the latitudinal univariate fatty acid comparisons, the significance level of the tests was adjusted using a Bonferroni-correction for multiple tests (where α = 0.05/10 tests = 0.005). Univariate fatty acid analyses were run in R version 4.1.0 (R Development Core Team 2021).

Probability of blue coloration: all surveys

Lingcod were caught across a similar range of depths for all regions except Washington, where fishing locations were deeper on average than other regions (Figure S2). To avoid the confounding effects of region and depth, we restricted analyses to fish caught at depths shallower than 125 m, resulting in removal of an additional 223 fish. To test hypotheses regarding potential drivers of blue flesh, we performed separate logistic regression analyses on the prepared data (n = 1798 lingcod; hereafter ‘full dataset’) and a subset of the full dataset for which fatty acid concentrations were measured (n = 163 lingcod; hereafter ‘FA only dataset’). For the full dataset, we used logistic regression to estimate the probability of blue coloration as a function of geographic region, depth, length, and sex. Region and depth were included in the models to evaluate whether latitudinal and depth distribution differences among individuals accounted for variation in coloration. Sex and length were included as proxies for sex- and size-based differences in diet and physiological status (e.g., reproductive hormones, growth hormones, and demographic shifts in diets) that are potential determinants of blue flesh. For the FA only dataset, five full models were constructed, each including region, depth, length, sex, and one of the following fatty acid variables: total fatty acid concentration, ω-3:ω-6 ratio, 18:1ω-9, 16:1ω-7, or 18:1ω-7 (concentrations). Total fatty acid concentration was included because of the significant differences found between blue and brown fish. The ‘ω-3:ω-6’ factor was included because of hypotheses generated from the independent southern Oregon dataset and because ω-6 fatty acids, particularly ARA, may be an important indicator of nutritional status in fishes (Sargent et al. 1999; Bell and Sargent 2003). The other three biomarkers were included because they were identified from the SIMPER analysis as the most influential fatty acids explaining differences between blue and brown fish (i.e., each contributing > 10%; see “Results”). Because of the relatively low sample size in the FA only dataset, including all potential fatty acid predictors in the full model would have resulted in a very large number of candidate models and can result in overfitting the data without contributing meaningful biological inference (Burnham and Anderson 2002).

For each of the datasets, multiple logistic regression models were fitted to the binary response of blue (1) or brown (0), with a logit link function using the ‘glm’ function in the R package ‘MASS’ (R Development Core Team 2021). Since the units of measurement and scales of predictor variables varied substantially, we standardized the factors prior to fitting the regression models (Gelman 2008) using the ‘rescale’ function in the R package ‘arm’ (R Development Core Team 2021). Rescaling predictor variables also allows for direct comparison of estimated coefficients (Gelman 2008). The sets of candidate models included the full models and all possible combinations of potential predictors, as described above for each dataset. Model selection was performed using Akaike’s Information Criterion corrected for small sample size (AICc), where models with AICc values within two units of the lowest AICc across all models were included in the set of top models (Burnham and Anderson 2002). We evaluated the relative importance of explanatory factors using parameter weights, calculated as the sum of Akaike model weights across all models including a given parameter (Burnham and Anderson 2002). Parameter weights range from 0 to 1, where values approaching 1 signify a strong weight of evidence that the predictor should be included in the best-fit model (Burnham and Anderson 2002).

Results

After removing fish with no length measurement, unidentified sex or color, and colors with mismatched external and internal color designations, we were left with 2,021 fish for analysis (Table S1) in the latitudinal survey dataset. Of these, 286 (14%) were blue in color. Most blue fish were female (81%) and the proportion of blue lingcod within each region ranged from 3% (Washington coast) to 21% (northern California) (Fig. 2).

Southern Oregon survey

The large majority (90%) of the lingcod in the initial southern Oregon survey with blue flesh were female, differing strongly from an even sex ratio (Pearson’s Chi-squared test; χ2 = 20.18, df = 1, p < 0.0001). Of females, 55% were blue-fleshed and 45% were brown, whereas only 5% of males were blue. Multivariate fatty acid composition in the southern Oregon dataset did not differ by color, sex, or their interaction (2-way PERMANOVA; color: Pseudo F = 0.83, df = 1, p = 0.41; sex: Pseudo F = 2.09, df = 1, p = 0.11; color x sex: Pseudo F = 0.96, df = 1, p = 0.36; no nMDS presented). The southern Oregon lingcod dataset was limited in the sample size of blue fish. Since almost all blue fish were female, preliminary analyses of univariate fatty acid summary categories focused on female fish. Blue female lingcod had a higher ω-3:ω-6 ratio (t test = 2.23 df = 26.70, p = 0.035), lower total log(ω-6) fatty acids (t test = -2.25, df = 21.93, p = 0.034), and lower log(20:4ω-6) (ARA) content (t test = 2.11, df = 21.52, p = 0.046), but these initial comparisons were not statistically significant when Bonferroni-correcting for multiple tests in the Oregon dataset. These preliminary findings with small sample sizes were used as a pilot study to inform the later analysis and are, therefore, interpreted with caution.

Latitudinal survey fatty acids

Lingcod multivariate fatty acid concentrations differed significantly between regions and fish color, but not sex; there were no significant interactions among any terms and the full model explained ~ 26% of the total variation (3-way PERMANOVA; Table S3). Multivariate fatty acid dispersion did not differ among regions (PERMDISP; F6, 166 = 1.39, P (perm) = 0.57) or sexes (PERMDISP; F1, 171 = 0.93, P (perm) = 0.45), but blue and brown fish differed in their dispersion (PERMDISP; F1172 = 9.42, P (perm) = 0.012). The nMDS ordination (Fig. 3) showed generally lower multivariate fatty acid dispersion within blue-fleshed lingcod, with tighter grouping on the left side of the plot (except for 3 outliers, which were included in the PERMDISP test).

Fig. 3
figure 3

Non-metric multidimensional scaling (nMDS) plot showing multivariate dispersion of blue- and brown-fleshed lingcod. The analysis includes 173 fatty acid samples from the West Coast latitudinal dataset. Data are Log(X + 1) transformed fatty acid concentrations (mg FA g−1 dry weight)

SIMPER identified that 18:1ω-9, 16:1ω-7, 18:1ω-7, 20:5ω-3, 22:6ω-3, 16:0, 14:0, 20:1ω-9, 20:4ω-6, 18:0, and 22:5ω-3, in that order, were the most important fatty acids for differentiating blue and brown lingcod (see SIMPER results, Table 1). Blue fish (pooled across sexes and regions) had significantly lower concentrations of log 16:1ω-7 (t test = 4 .09, df = 169.24, p = 0.0001), log 18:1ω-9 (t test = 4 .63, df = 171, p < 0.0001), log 18:1ω-7 (t test =  4.35, df = 168.72, p < 0.0001), log 20:4ω-6 (t test =  4.05, df = 166.51, p = 0.0001), log 20:5ω-3 (t test =  4.37, df = 164.39, p < 0.0001), log total ω-6 (t test =  4.67, df = 159.81, p < 0.0001), and log total fatty acids (t test =  3.47, df = 170.37, p = 0.0007) (Fig. 4). There were no differences in the concentrations of log 22:6ω-3 (t test = 0.42, df = 166.99, p = 0.67) or log total lipids (t test =  1.57, df = 160.37, p = 0.12; not pictured). Blue fish had a significantly higher ω-3:ω-6 ratio (t test = 4.31, df = 170.18, p < 0.0001) than brown fish (Fig. 4).

Table 1 Multivariate similarity percentages (SIMPER) analysis on the log(X + 1) transformed fatty acid concentration data (mg FA g−1 dry wt.) was used to identify the fatty acids that were most influential in separating the blue and brown fish visualized in the nMDS (Fig. 3)
Fig. 4
figure 4

Data (colored points) are jittered within each category to reduce overplotting. Fatty acids were included based on SIMPER results and a priori hypotheses. ac Key mono-unsaturated fatty acids; df key polyunsaturated essential fatty acids; g the sum of all ω-6 fatty acids; h the ω-3:ω-6 fatty acid ratio; and i total fatty acids. All comparisons differ significantly except for panel f (indicated by “*”; Bonferroni-corrected t tests, where α < 0.005)

Probability of blue coloration

We modeled the probability of blue flesh as a function of region, capture depth, total length, and sex using the full dataset. The best model included all of these variables (Table 2), which all had high parameter weights (> 0.8). Sex had the greatest effect on probability of blue coloration, followed by length, depth, and region (Table S4). Based on the signs of estimated coefficients from the best-fit model, lingcod that were female, relatively smaller, and relatively shallower in bottom depth at capture had a higher probability of blue coloration (Fig. 5). There was no apparent latitudinal trend in the proportion of lingcod showing blue coloration and regional differences were driven largely by the Washington coast, where the predicted probability of blue coloration was smallest, and fewer samples were collected in shallower waters near the coast.

Table 2 Model diagnostics from analyses of the full dataset and the fatty acid ‘(FA) only’ dataset
Fig. 5
figure 5

Distributions of lingcod length across 25 m depth bins from the full West Coast latitudinal dataset. Separate boxplots are shown for each combination of color and sex. Boxes show the interquartile range (IQR), with the median shown as a line within each box. Whiskers show ± 1.5*IQR. Outliers were removed from the figure for ease of visualization

Fatty acid concentrations were relatively unimportant in determining the probability of blue coloration (Table 2). Analysis of the fatty acid-only dataset revealed that eight models performed equivalently (∆AIC ≤ 2) and, when considered together, included all potential predictors. However, consistent with the analysis of the full dataset, the most parsimonious model (∆AIC = 0) included only sex, length, and depth (Table 2). Sex and length were the most influential predictors (parameter weights = 1) with the greatest effect sizes, and depth was a relatively important predictor with parameter weight = 0.78 (Table S4). As in the full model, lingcod that were female, relatively smaller, and residing in relatively shallower depths had a higher predicted probability of blue coloration compared to males, relatively larger individuals, and those residing in relatively deeper depths (Table S4). Region, total fatty acid concentration, ω-3:ω-6 (selected a priori), 18:1⍵-9, 16:1⍵-7, and 18:1⍵-7 (the top three fatty acids from the SIMPER) were relatively unimportant in explaining variation in occurrence of blue lingcod, with effect sizes not different from zero (Table 2) and parameter weights < 0.3 (Table S4).

Discussion

The question of why some lingcod are blue has long been of interest to seafood workers, fishers, scientists, and naturalists. More than an ecological curiosity, this phenomenon has commercial implications and possible evolutionary consequences that deserve further investigation. Our study was the first to use an extensive dataset from across the species range to improve our understanding of spatial and demographic patterns in the prevalence of blue fish and potential dietary and physiological drivers of blue flesh. We found that the probability of having blue flesh was highest for fish that were female, caught in shallower water, and smaller in total length. Sex is a previously overlooked factor for blue color in lingcod, given that > 80% of blue fish across this large geographic range were female. Incidence of blueness was 14% coastwide and varied by region (3–21% of all fish) but was also confounded by sex ratios of fish caught among regions. In lingcod, smaller females are typically found inshore and migrate to deeper waters at larger sizes, returning only for reproduction (Low and Beamish 1978); however, some reproductively mature female lingcod remain in shallow nearshore habitats. For example, an acoustic telemetry study in central California found that small females had the highest residency within the nearshore array (Greenley et al. 2016). Considering our results showing that smaller, shallower, female fish are more likely to be blue, it is possible that remaining in nearshore environments is associated with a physiological response or a dietary difference that is less common in larger, deep dwelling females and cohabitating males. Recent research has shown that male blue lingcod also carry heavier burdens of metazoan parasites than do brown males, brown females, or blue females (Wood et al. in revision), and that those blue males are in poorer condition (low Fulton’s K condition factor) than other sex-color combinations. It may be that starvation stress causes the release of compounds causing blueness and also weakens immune defenses, resulting in higher parasite burdens for males only because males are already immunocompromised by male sex hormones (Folstad and Karter 1992; Roberts et al. 2004; Foo et al. 2017).

Biliverdin has long been identified as the underlying cause of blue coloration of serum and other body parts in fishes (Fang and Bada 1990). We observed that lingcod with blue flesh reliably had visibly blue-colored serum. Biliverdin is found in the bile and serum of a diverse lineage of animals. In mammals, it is rarely detectable, except during hepatic dysfunction, as it is enzymatically converted to yellow/orange bilirubin; however, it is the primary bile pigment in many avian, reptile, amphibian, and fish species (Cornelius 1984). In fishes, biliverdin may serve a role in external coloration (Gagnon 2006), starvation, or as a signal of stress (Love 1980). For example, attempts to replace traditional fish meal with plant-based proteins in aquaculture of red sea bream (Pagrus major) and yellowtail (Seriola quinqueradiata) have resulted in a taurine deficiency, linked to increased levels of biliverdin in hepatic fluids and resulting in a condition known as green liver (Takagi et al. 2010). Starvation has also been shown to lead to increased levels of serum biliverdin in a number of fish species (Fang 1987). For example, in Atlantic cod (Gadus morhua) it has been shown that with experimental starvation there is an excess production of blue-green bile in the gall bladder, which leaks into the intestines (Love 1980). In our initial exploratory analyses of the lingcod data, we did not identify total lipid content or the hepatosomatic index (HSI), which are common indices of condition or starvation, as important factors in predicting blueness; however, using a different subset of samples, Wood et al. (in revision) report that blue lingcod of both sexes have lower HSI than brown individuals. Other hypotheses for why biliverdin accumulates in the tissues have been advanced, but not tested. For example, Yamaguchi et al. (1976) posited that the biliverdin in blue-fleshed sculpins (Pseudoblennius percoides) might be stored in the flesh before being transported to skin and eggs for the purpose of camouflage, ultra-violet protection, or through involvement with reproductive activities, when needed. We collected additional samples of lingcod and their serum in Puget Sound, Washington, (n = 70) during opportunistic research efforts during spring–summer of 2018 and 2019 in an attempt to analyze biliverdin and serum coloration. Unfortunately, we were not able to replicate other fish serum biliverdin analysis methods, despite the very clear blue coloration of the serum of blue-fleshed fish (Figs. S3, S4). The serum pilot analysis is described in further detail in the Supplemental Materials. Further research is needed to characterize conditions under which biliverdin production and excretion is modulated in lingcod.

The association between fatty acid concentrations and lingcod blueness was nuanced. Multivariate fatty acid signatures differed between blue and brown lingcod, but an nMDS visualization showed considerable variability and logistic regressions suggested relatively little support for individual fatty acids as predictors of blueness. SIMPER analyses identified that the most important fatty acids for differentiating blue and brown fish were 18:1ω-9, 16:1ω-7, 18:1ω-7, 20:5ω-3 (EPA), 22:6ω-3 (DHA), and 16:0, which accounted for 73% of the cumulative variation in the SIMPER analysis. Most of these fatty acids had small absolute differences in terms of total mg fatty acids g−1 of tissue, and while most of these differences were statistically significant, without additional experimental work on lingcod physiology, it is unknown if these differences are biologically significant. The five fatty acid variables (18:1ω-9, 16:1ω-7, 18:1ω-7, ω-3:ω-6, and total FA) tested in logistic regression models were relatively unimportant in explaining the probability of blue flesh, when effects of sex, length, and depth were accounted for. Several of the most important fatty acids for differentiating blue and brown fish (16:0, 18:1ω-9, and 20:1ω-9) have also been described as the primary sources of metabolic energy in fishes (Tocher 2003). Concentrations of each of these fatty acids were lower in blue lingcod, which might indicate that blue fish are obtaining relatively more energy from catabolization of energy stores. The mono-unsaturated fatty acid 18:1ω-9 (oleic acid) was the most important fatty acid in the comparison of blue and brown lingcod, based on both the multivariate and regression analyses. Oleic acid is not a particularly valuable trophic biomarker of specific sources because it is common in most marine organisms, and can commonly be synthesized de novo (Sargent et al. 2002; Tocher 2003). However, high levels of oleic acid or the ratio of 18:1ω-9 to 18:1ω-7 have been identified as possible markers for carnivory in some zooplankton (Brett 2009); it is unclear if this relationship translates to fishes. Evaluating all 11 fatty acids identified with SIMPER showed that the concentrations of all 4 mono-unsaturated and 3 saturated fatty acids were lower in blue fish. DHA and EPA, which are considered physiologically essential fatty acids in marine fishes (Tocher 2003), were the most important PUFAs in separating blue and brown lingcod, despite the fact that DHA did not differ significantly between blue and brown fish.

Total fatty acid content was significantly lower in blue fish (Fig. 4i), which have 17.5% lower mean total fatty acid concentrations (Table S2) compared to brown fish. Total fatty acids are sometimes used as an indicator of fish condition (e.g., Pethybridge et al. 2015); whether these observed differences are physiologically significant in lingcod is unknown. In addition, we found that blue-fleshed fish had significantly lower concentrations of ω-6 fatty acids, which are precursors for the synthesis of eicosanoids, which modulate inflammation and facilitate the immune response in fishes (Bell and Sargent 2003; Arts and Kohler 2009). Excess eicosanoid synthesis in fishes is known to be initiated by stress or other stimuli, triggering the oxygenation of these fatty acids through various pathways, resulting in ARA-derived eicosanoids which promote inflammation and EPA-derived anti-inflammatory eicosanoids (Arts and Kohler 2009). Blue fish had lower content of 20:4ω-6 (ARA) and EPA, which are the primary fatty acids from which eicosanoids are derived (Sargent et al. 1999; Arts and Kohler 2009). Our observation of lower ω-6 content in blue fish is still just an emergent coincident pattern; it is unclear if lower ω-6 content is causatively linked to the blue flesh, particularly given the weak effects of fatty acids after accounting for size, sex, and depth.

The fatty acids of marine consumers are driven by a combination of diet, physiology, and environmental conditions (Dalsgaard et al. 2003; Xu et al. 2020). In this study, multivariate fatty acid composition differed among fish across regions and flesh colors, but not between sexes. Sex was the most important factor for predicting the probability of blue flesh in the logistic regression analysis, possibly swamping out effects of fatty acids. Further investigation of the role of sex steroids and other factors, such as parasite load, in explaining phenotypic variation in fish coloration may help to better understand why the prevalence of blue flesh varies among sexes. The multivariate results did not differ when we evaluated the fatty acids of only female lingcod. The lack of differences in fatty acids between males and females suggests that lingcod of both sexes have generally similar diets, or that both sexes experience the same physiological response to the underlying driver of blueness. The reason for the three blue outliers on the right side of the nMDS ordination is unknown. Data exploration did not identify any particular explanatory factor (fish length, depth captured, region, sex, etc.) that could account for their anomalous position relative to other blue fish. Despite this variation, the PERMDISP analysis still indicated that brown fish were more variable than blue fish in their fatty acid composition. The fatty acids of lingcod differed across the large geographic range we investigated, which may indicate differences in trophic partitioning of lingcod in the different major eco-regions of the Eastern Pacific. Region did not interact with color in the fatty acid analysis, indicating that blue fish were depleted in key fatty acids compared to brown fish in all regions sampled. Future work using additional tools for trophic inference, including stomach contents and stable isotope analysis, should leverage these initial findings to see if there are multiple lines of evidence to support the fatty acid results, which suggest divergent trophic ecology of the fish across regions. Our analysis was based on the commonly used total lipid extraction on muscle tissue, but analysis of only neutral lipids may be preferable for assessment of diet in a predator (Dalsgaard et al. 2003). We used muscle tissue, as this was consistently collected from the same part of the fish from all lingcod in the latitudinal study; but future work could also look at other tissues (e.g., adipose tissue or liver) which might better reflect recent diets (Budge et al. 2006).

In addition to the scientific interest in the physiological, ecological, and evolutionary roles that heterogeneous coloration patterns play in animal populations, this topic is of profound interest in the general public. Albino individuals among fish and wildlife species garner media stories in the popular press (Doughton 2007; Sanchez 2019; Dickson 2020). For lingcod, blue individuals are the subject of many forum conversations, and there are many opinions on what causes this color variation among anglers (e.g., http://www.theoutdoorline.com/forum/; https://www.sportfishingbc.com/forum/index.php). We recently received a confirmed report of a very rare bright-pink lingcod color morph caught by a recreational angler out of Newport Oregon in 2020, which had pink external and internal flesh coloration (unpublished observation reported to AWEG by Oregon Department of Fish and Wildlife). In the marketplace, the coloration of fish fillets also plays a role in the price and demand of seafood products, with blue lingcod fillets “always selling faster than white counterparts” (Bland 2014). This interest has the potential to lend itself to community (‘citizen’) science opportunities. Integrating data collected by recreational anglers on the sex, size, color, depth of capture, and stomach contents for this popular sportfish would not only provide new data across unsampled areas, but would also provide a framework for collaboration between science and stakeholders, which could improve trust in the management decision-making process (Enquist et al. 2017; Venturelli et al. 2017; Andrews et al. 2019).

Our analysis has provided new insights into the prevalence and factors associated with blue flesh for a very large sample of lingcod, across their geographic range. This is, to our knowledge, the largest such analysis of fish flesh coloration ever published. Most published work on color variants in other fishes highlight relatively rare or novel occurrences. These initial results draw attention to several outstanding mysteries. Why is it that female fish are more likely to have blue flesh? Our analysis identifies that sex has been an otherwise overlooked factor that is relevant for blue color in lingcod, given our discovery that > 80% of blue fish are also female. Does reproductive status, which may affect males and females differently, affect the blue coloration? A common hypothesis among fishers and biologists that blue-colored flesh in lingcod is due to diet (Love 2011). Unfortunately, our findings do not resolve this mystery. The differences in fatty acids among blue and brown fish can be taken as support for this hypothesis; however, it is still possible differences in fatty acids are a reflection of some other physiological process besides diet, and this question will not be resolved without lab experiments.