Way out there: pathogens, health, and condition of overwintering salmon in the Gulf of Alaska

Samples of Pacific salmon were collected during the 2019 International Year of the Salmon GoA expedition in February and March 2019 onboard the Russian research trawler Prof. Kaganovskiy. Sixty 1 h trawls accompanied by oceanographic sampling were performed along a grid of stations separated by 1 degree of latitude or 1.5 degrees longitude (approximately 110 km apart), and 422 salmon were captured over the course of the expedition (Supplementary Material, Fig. 1, Supplementary Material, Table 1). Subsampled salmon from all species were dissected in a clean environment within one hour of capture (Supplementary Material, Fig. 2). Notes on gross pathologies were collected during dissections. Presence of nematodes in organs or the peritoneum was noted on a nonspecies-specific level; no other macroscopic parasites were observed. Tissue samples were preserved in RNAlater (Thermo Fisher Scientific, MA, USA) for nucleic acid extraction as well as in 10% neutral buffered formalin for histopathology.

Genetic stock identification for coho and sockeye salmon was performed by the Department of Fisheries and Oceans Canada, Pacific Biological Station Molecular Genetics Laboratory as described by Beacham et al. (2010, 2020).

Oceanographic data were collected at each station with a 24-position rosette equipped with CTD sensors as described by Pakhomov et al. (2019). Turbidity, fluorescence, and oxygen saturation were measured, and water samples were collected for assessing salinity, chlorophyll, and macronutrients. To survey zooplankton communities two Juday nets as well as one Bongo net were deployed as described by Pakhomov et al. (2019).

Calorimetric data on the energy content of salmon individuals in this study was provided by the National Oceanic and Atmospheric Administration, Alaska Fisheries Science Center, Auke Bay Laboratories in Juneau, Alaska, USA, for 46 Sockeye. In brief, tissue samples were weighed, dried with their skin on at 135 °C, homogenized, and analyzed using bomb calorimetry (kJ/g dry mass) (Siddon et al. 2013).

Tissue samples from gill, heart, kidney, and spleen were homogenized using TRI-reagent (Ambion Inc., Austin, Texas), and 1-bromo-3-chloropropane was added to the homogenate. Total RNA was extracted from the aqueous phase using the Total RNA Isolation kits (Ambion Inc., Austin, Texas) on a Biomek FXP liquid handling instrument (Beckman-Coulter, Mississauga, Ontario, Canada) (Miller et al. 2011; Jeffries et al. 2014). RNA quality was assessed after DNase treatment by spectrophotometry and RNA was normalized to 62.5 ng/μL so that 1 μg of RNA was used for cDNA synthesis (SuperScript VILO MasterMix, Life Technologies).

DNA was extracted from the organic/interphase of TRI-reagent using a high salt TNES-urea buffer (Asahida et al. 1996) followed by the BioSprint 96 DNA extraction kit (Qiagen, MD). DNA quantity and quality were assessed by spectrophotometry on a Beckman Coulter DTX 880 Multimode Detector (Brea, CA, USA). Samples were normalized to 62.5 ng/μL.

For infectious agent monitoring, cDNA from all pooled organs was mixed with equal amounts of pooled DNA extract from all organs to 1.25 μL final volume. Samples were pre-amplified with primer pairs for all 48 infectious agent assays (Table 1) in a 5 μL reaction using TaqMan PreAmp Master Mix (Life Technologies), following the BioMark protocol, to increase sensitivity of the small assay volume (7 μL) on the dynamic arrays. Unincorporated primers were eliminated using ExoSAP-IT PCR Product Cleanup (MJS BioLynx Inc., Ontario, Canada). Samples were diluted 1:5 in DNA Suspension Buffer (TEKnova, Hollister, California). An assay mix containing 9 μM artificial positive control clones (labelled with fluorescent dye VIC) allowed for the detection of contamination. For host gene expression monitoring, an equivalent procedure was performed on cDNA from gill tissues only, targeting 89 host genes individually. A serial dilution of pooled gill cDNA was used to assess assay efficiency across runs.

The qPCR assays and individual samples were loaded onto 96.96 dynamic arrays (Fluidigm, San Francisco, CA, USA) and run on the BioMark^TM HD platform. The same distribution of assays was used for each array and samples from different dates and locations were stratified among arrays. The Fluidigm 2 × Assay Loading Reagent was mixed with primer pairs and probes to 5 μL per well. A 5 μL sample loading mix was prepared using 2 × TaqMan Gene Expression Master Mix (Life Technologies), 20 × Gene Expression Sample Loading Reagent, and 2.7 μL of pre-amplified cDNA. The reaction mixes were added to the assay and sample inlets of the dynamic array as per the manufacturer’s protocols and loaded into the chip by an integrated fluidic circuit controller HX (Fluidigm). PCR was performed under the following conditions: 50 °C for 2 min, 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min.

In addition to known pathogens, we incorporated assays for newly discovered putative viruses of unknown pathogenic potential. These viruses originate from an unpublished polyA amplified metatranscriptomic sequencing library from 20 Chinook salmon targeting unknown infectious agents. We screened these libraries using a translated blast search (see Mordecai et al. 2019 for methods) and found short contigs (Genbank accessions MW373508-MW373514, Supplementary Material, Table 2). These contigs showed protein homologies to hanta-like, rhabdo-like, picorna-like, and qin-like viruses, as well as contigs with high sequence homology to an unpublished viral contig SalmovirusWFRC1 (NC_034441). Additionally, we included an assay for a sequence variant of Pacific salmon Nidovirus (PsNV; Mordecai et al. 2020) and a co-infected bafini-like virus (Supplementary Material, Table 2).

Cycle threshold (CT) for each assay was determined using the BioMark Real-Time PCR analysis software (Fluidigm). For each infectious agent assay, samples with detection in only one duplicate were treated as negatives and duplicate values were averaged. Samples contaminated by high load controls (indicated by VIC positives) were removed. Amplification curves of all assays were visually assessed for irregularities and consistency between replicates. R statistical software (R Core Team 2017) was used to calculate the efficiencies for each assay using the slope of a regression between CT values and serial dilutions of the artificial positive control standards. We removed values that were not within the linear relationship, often either the lowest or highest RNA concentrations, to improve accuracy of assay efficiency estimates and r² values. Only assays with an amplification factor of 1.80–2.20, an r² value of ≥0.98, and with typical shaped amplification curves were used in analyses (Larionov et al. 2005). Minimum averaged CT values indicating infectious agent detection with high statistical certainty for each specific infectious agent assay (95% confidence limit of detection), were defined by Miller et al. (2016). Infectious agent prevalence was calculated as the percentage of individuals testing positive for a given infectious agent. All infectious agents found within one host were summarized as a single variable termed “Relative Infection Burden” that takes into account the infection load of all detected pathogens in an individual compared with the population average (Bass et al. 2019). Infectious agent load is the number of copies of a given infectious agent in an individual testing positive. To examine differences between high-sea samples and coastal populations we compared the GoA data with baseline data from coastal British Columbia based on 11 790 wild or nonhatchery marked salmon of all species and age classes sampled between 2014 and 2019 between the Juan de Fuca Strait in the south and waters at the Alaskan border near Dixon Entrance. We calculated location (coast vs GoA) and species-specific relative infection burden (RIB) and Shannon Weaver diversity of infectious agents and compared prevalence of specific agents on a species-specific level using Fisher’s exact test (Clarke and Warwick 2001).

To determine the primary stressors experienced by salmon in the GoA, we deployed salmon Fit-Chips that utilize curated panels of 89 host genes (biomarkers) to detect transcriptional responses to stressors in gill tissue on the same nanofluidics qPCR platform described above (Miller et al. 2011, 2017; Akbarzadeh et al. 2018; Houde, Akbarzadeh, et al. 2019; Houde, Günther, et al. 2019). Physiological states are recognized based on co-expression of curated biomarker panels that have consistently segregated stress and disease states in challenge studies. For the GoA samples, we applied biomarker panels for hypoxic stress, thermal stress, osmotic stress, general stress, and viral disease development (genes expressed in response to active viral infection), as well as imminent mortality (over-expressed in salmon experiencing mortality within 72 h) and mortality-related (associated with poor long-term survival) markers (Miller et al. 2011, 2017; Akbarzadeh et al. 2018; Houde et al. 2019; Houde et al. 2019). We also included biomarkers associated with different branches of immune stimulation (over-expressed in diseased individuals with known pathogens) and with inflammation (individuals showing pathological signs of inflammation; Table 2). All biomarkers have been assessed for efficiency of amplification across all salmon species, but development of the panels used Chinook and sockeye. Applications across four salmonid species herein offers our first glimpse into recurring patterns of stress- and disease-related gene expression patterns across species co-inhabiting offshore waters of the North Pacific.

Host genes assays were run singularly on cDNA from gill tissues and included three reference genes for normalization (Miller et al. 2016; Teffer et al. 2017). Host gene assay efficiencies were calculated using the serial dilution of pooled pre-amplified host cDNA run on each dynamic array (Miller et al. 2016). Expression heatmaps were visually assessed for failed assays or samples, samples with low expression of reference genes were removed, and failed assays were assigned the mean of the respective species. Samples with less than 55 ng/μL cDNA were excluded from the analysis. Salmon gene CTs were normalized between runs one species at a time using calibrator samples, converted to relative expression by normalizing against the average of the best two out of three reference genes as determined by normfinder, and the relative fold gene expression was calculated using the ddCT method (Livak and Schmittgen 2001; Jensen and Ørntoft 2004). Assays that failed in more than 50% of individuals for the respective species were excluded from analysis.

To gain an overview of gene expression and cluster individuals into groups, expression profiles were visualized as heatmaps using the package ComplexHeatmap in R (Gu et al. 2016). Heatmaps were augmented with pathogen and significantly co-varied metadata between gene expression clusters as determined by analysis of variance (ANOVA) and t-test analysis in R (base, stats). To determine the dominant stressors experienced at a population level, we compared the expression of all genes in all individuals of a given species using principal component analysis (PCA). For multidimensional data, PCA identifies the dominant axes (or dimensions) of variation, allowing quantitative interpretation of differential expression among individuals through “ordination” at reduced dimensionality. We deployed the prcomp function in R (base, stats). For visualization, we depict all individuals in the first four dimensions of the PCA, as well as showing the top 20% of genes responsible for ordination in the depicted dimensions as determined by the ordiselect function of the package goeveg in R (cran.r-project.org/web/packages/goeveg/). Three outlier individuals that had dissection comments suggesting severe damage during capture explaining the aberrant gene expression profiles were excluded from the analysis. To interrogate correlations of gene expression with infectious agent, physical, and oceanographic data at site of capture, this information was ordinated onto the PCA plots using the envfit function of the vegan package in R (cran.r-project.org/web/packages/vegan/index.html; Table 3). This function scores correlation of data with the given ordination dimensions and provides a quantitative directional vector depicting this correlation. For visualization, only significant vectors with p < 0.05 after 999 permutations are displayed with the metadata name indicating the tip of the scaled arrow segments (Supplementary Figures). We summarized the vectors of all genes belonging to a specific biomarker, and we evaluated the correlation of this superimposed data with the genes driving the ordination using ordiselect in R and the top 20% of genes showing significant correlation in expression with the data (see figures in Results section).

In addition to PCA ordination of gene expression, we also deployed nonmetric multidimensional scaling (NMDS) to describe the different pathogen profiles carried by individuals using the metaMDS function of the vegan package in R. This ordination approach is preferable to PCA for the pathogen data, as it produces an ordination based on abundance rank order, rather than absolute values, which is better able to deal with missing data (in this case absence of pathogen detections). To find the correlations of pathogen profiles with gene expression, physical, and oceanographic data, we used the same approach as described for PCA data above (Table 3).

To investigate stressors of salmon in the GoA, we compared the expression of all genes from all biomarker panels across all individuals of the same species. First, we visualized gene expression using heatmaps, also displaying pathogen detections as well as co-varying metadata (Fig. 3). Hierarchical clustering of gene expression allowed us to identify clusters of salmon showing similar expression patterns (Fig. 3). In chum, clusters four and five showed markedly reduced overall gene expression that was associated with elevated RIB and lower biomass of hydromedusa at capture location, the primary prey of chum salmon (Somov et al. 2019), as well as lower levels of dissolved oxygen (Fig. 3; Supplementary Material, Fig. 5). Further, temperature at site of capture was also significantly associated with overall gene expression, with warmer temperatures correlating to higher gene expression (Fig. 3; Supplementary Material, Fig. 5). Similarly, in sockeye, elevated temperature, and prey availability (e.g., small zooplankton) was associated with a global increase in gene expression (Fig. 3; Supplementary Material, Fig. 6). Condition factor K was significantly covaried between clusters in sockeye (Fig. 3; Supplementary Material, Fig. 6). Coho salmon showed no large-scale changes in gene expression and clusters differed primarily in the response of individual biomarkers to RIB and prey availability (Fig. 3; Supplementary Material, Fig. 7). Pink salmon also showed large-scale changes to gene expression associated with RIB, prey availability, and temperature, but interestingly high values of these factors were associated with reduced global gene expression rather than an increase as had been seen in chum and sockeye salmon (Fig. 3; Supplementary Material, Fig. 8).

Next, we performed a PCA to ordinate gene expression profiles of individual salmon. We focused on the first four principal components to identify dominant biomarker panels driving differential gene expression. We then tested observational data on salmon health and condition as well as oceanographic data for correlations with principal components and plotted significantly correlated data scaled and directional on the ordination plots to depict the direction of correlation. In the last step, we queried what genes showed changes in expression correlated with the superimposed data by using a Euclidean distance-based approach or plotted a vector summarizing of all genes of a biomarker panel, respectively (Fig. 4; Supplementary Material, Fig. 9). By visualizing the hierarchical clusters identified earlier, this allowed us to identify the environmental factors and pathogens associated with the differential gene expression providing a population scale overview of stressors of overwintering salmon (Fig. 4; Supplementary Material, Fig. 9).

Differential gene expression in chum salmon was primarily driven by variations in biomarkers for inflammation (MMP13, NAPEPLD2, TXN, GILT), immune stimulation (SAA, CD83, IFNa), mortality-related (C7, P_RAS), viral disease development (VDD) biomarker panels (HERC6, IFIT5, IFI44, VAR1), followed by imminent mortality and hypoxia (TAGLN3, CDKN1B; Fig. 4a; Supplementary Material, Fig. 9a). Along PC1, these factors explained 36% of the variation in gene expression. Chum clusters four and five showed lower gene expression across all biomarker panels and clustered on the positive end of PC1 (Figs. 3 and 4a; Supplementary Material, Fig. 9a). Inflammation (MMP13, NAPEPLD2, TXN) and immune stimulation biomarkers (SAA) contributed the negatively to PC2 (10.8% explanatory power), while hypoxia and imminent mortality biomarkers (CDKN1B) contributed positively (Figs. 3 and 4a; Supplementary Material, Fig. 9a). RIB as well as nematode prevalence was correlated with lower overall gene expression in individuals of cluster four and five, but positively associated with inflammation (MMP13, NAPEPLD2), immune stimulation (SAA), and VDD (HERC6, IFIT5) markers on PC2. Conversely, pathogens P. pseudobranchicola (pa_pse) and S. destruens (sp_des) were negatively associated with these immune response markers along PC2 (Fig. 4a; Supplementary Material, Fig. 9a). Biomass of hydromedusae (Medu) and other the prey of chum (small zooplankton: Zoo_S) was positively correlated with global upregulated gene expression along PC1 and lower expression of the immune response markers associated with PC2 (Somov et al. 2019), while being directly opposed to RIB and nematode prevalence across PC1 and PC2 (Fig. 4a; Supplementary Material, Fig. 9a). Principal components three and four (explaining 10% and 6.2%, respectively) were driven by the same genes driving PC1 and PC2; however, inflammation and VDD markers were opposing each other along PC3, with inflammation driven by individuals of cluster four, that showed enlarged gallbladders (a sign of prolonged low stomach fullness) and larger size (Mass) and smaller individuals at higher temperature associated with VDD expression (Fig. 4b; Supplementary Material, Fig. 9b). PC4 was driven by opposing trends of immune stimulation and hypoxia biomarkers, primarily associated with zooplankton (euphausiids and medium size zooplankton; Figs. 3 and 4b; Supplementary Material, Fig. 9a).

Sockeye showed similar patterns to chum salmon with two clusters (one and four) showing reduced overall gene expression associated with the positive end of PC1 (43.8%). The primary drivers associated with these global expression changes were the biomarker panels immune stimulation (B2M, HEP, IGMs, CD83, SAA), inflammation (IL_17D, ES1), mortality related (SCG2, RPL6), VDD (HERC6, DEXH, MX, IFI), and a group of hypoxia genes (RRl, CLASPIN, KIF15, COX6B, RRM2) (Fig. 4c; Supplementary Material, Fig. 9c). These hypoxia genes were also major contributors to PC2 (13.8%) opposed by the general stress marker JUN_F3 (Fig. 4c; Supplementary Material, Fig. 9c). Globally lowered gene expression in sockeye clusters one and four was associated with lower abundance of small zooplankton (Zoo_S), pteropods (Ptero), and hydromedusae (Medu) along PC1, and to a lesser degree lower temperature at site of capture (Fig. 4c; Supplementary Material, Fig. 9c). Euphausiids (Euphaus) that were identified as the primary prey of sockeye, were correlated with the positive end of PC2, opposed to temperature, and showed increased expression of inflammation and immune stimulation markers, but showed weaker association with gene expression than other prey groups (Fig. 4c; Supplementary Material, Fig. 9c). The prevalence of the gill parasite Loma spp. (lo_sal) was associated with expression of inflammation and immune stimulation biomarkers along PC1 and PC2. Principal component three (7.7% exploratory power) saw a strong correlation of immune stimulation (SAA, IFNa, IGM) and inflammation biomarkers (IL_17D, MMP24, MMP13) with the parasites I. hoferi and P. pseudobranchicola, whereas inflammation (ES1, EPD) and imminent mortality markers (TAGLN3, RGS21) were associated with nematode prevalence (Fig. 4d; Supplementary Material, Fig. 9d). Fish with higher caloric content and better condition factor (K) were also associated with lower temperatures at site of capture and lower prevalence of pathogens (ic_hof, pa_pse) (Fig. 4d; Supplementary Material, Fig. 9d).

Differential gene expression in coho salmon showed a nuanced response of biomarker panels along the first two principal components where inflammation (MMP13, IL_11, NAPEPLD2, IL_17D), general stress (JUN_F3), immune stimulation (IL_1b, HEP, SAA, IFNa), and VDD (HERC6, GAL3) associated positively with RIB and fish of cluster four on the positive end of PC1 (21.2%; Fig. 4e; Supplementary Material, Fig. 9e). RIB was inversely related to the biomass of pteropods (ptero) that were the preferred prey of coho salmon in GoA in 2019 (Somov et al. 2019), with fish from cluster four experiencing the lowest pteropod biomass (Fig. 4e; Supplementary Material, Fig. 9e). Hypoxia biomarkers (COX6B, RRM2, CDKN1B) were correlated with the prevalence of the gill parasite Loma sp. (lo_sal) along PC1 and PC2 (17.5%) (Fig. 4e; Supplementary Material, Fig. 9e), specifically amongst individuals of clusters two, three, and five. Principal components three and four (12.4% and 7.1%, respectively) showed a global increase in expression that was associated with the size of fish (Mass) of individuals in cluster four, as well as an increased expression of VDD biomarkers related to P. pseudobranchicola (pa_pse) load (Fig. 4f; Supplementary Material, Fig. 9f).

Pink salmon of cluster one showed reduced global gene expression compared to other clusters, grouped along the negative spectrum of PC1 (64.6 %) and were associated with higher temperatures, higher RIB, and higher biomass of prey species (Fig. 4g; Supplementary Material, Fig. 9g). On the positive spectrum of PC1, clusters two, three, and four were associated with increased expression imminent mortality markers (CDKN1B, CBEBP, GPX) but were differentiated along PC2 (10.2%) with expression of hypoxia (COX_6B, RRI, GPX), and inflammation biomarkers (NAPEPLD2, IL_17D) were associated with cluster four, while clusters two and three showed increase expression of VDD (TRIM, GAL3, MX, VAR1, IFI) and immune stimulation markers (SAA, IGMs, HEP, CD83) that were associated with increased RIB, number of infectious agents, as well as the prevalence of the parasites S. destruens (sp_des) and P. psudobranchiola (pa_pse) (Fig. 4g; Supplementary Material, Fig. 9g). Principal component three (explaining 5.6%) showed elevated expression of VDD (GAL3, Mx, IFI, HERC6), immune stimulation (SAA, IL_15), and general stress genes (HSP90) along the positive end of PC3, which was correlated with larger individuals (Mass; Fig. 4h; Supplementary Material, Fig. 9h). Higher biomass of Euphausiids (Euphaus) along PC4 (4.1%) correlated with VDD expression and inflammation (EPD).

To highlight overlying trends of pathogens and environmental factors such as prey biomass and temperature, we plotted the biomass of primary prey species in relation to ocean temperature and RIB across the first two principal components of gene expression (Fig. 5; Supplementary Material, Fig.10). Since global depression of immune response genes (immune stimulation, inflammation, and viral disease development) effectively equals immunosuppression, we created the inverse vector of gene expression of said biomarker panels to depict this suppressed immune status. Indeed, immunosuppression showed an inverse relationship with the biomass of the primary prey species as well as a direct correlation with RIB in all species (Fig. 5). In chum and pink salmon this trend dominated gene expression along PC1 (Fig. 5). Coho showed a strong inverse correlation between primary prey biomass and RIB, but immunosuppression was only weakly associated with them along PC2, suggesting that large-scale changes in gene expression resulting in immunosuppression are subordinate to other factors relating to RIB (Fig. 5). In sockeye, gene expression patterns were more strongly associated with small zooplankton, rather than the primary stomach content which was Euphausiids (Fig. 5). Accordingly, lower biomass of small zooplankton was associated with immunosuppression and elevated RIB in sockeye along PC1 (Fig. 5). Coho and pink salmon that were primarily caught along the southern border of the distribution area and experienced the highest ocean temperatures and showed a strong correlation of immunosuppression and RIB with increased temperature (Fig. 5).

RIB of microparasites in the GoA was lower compared with coastal samples in all species except coho, which had a significantly higher relative infection burden due to high prevalence of the virus VER as well as the parasites I. hoferi and Loma sp. These two parasites and the bacterium Ca. B. cysticola were the highest prevalence pathogens in the GoA across all species.

Ichthyophonus hoferi was present at significantly higher prevalence and load in all salmon in the GoA compared with coastal areas. This common parasite causes systemic disease in marine fish and is thought to transmit trophically (Hershberger et al. 2002; Bass et al. 2017). This suggests that the GoA is a reservoir for this parasite and that piscivorous species acquire infection through their prey. Ichthyophonus hoferi detections in chum salmon, a species with low proportion of fish in its diet (1.8% in the study area) is surprising but suggests very high I. hoferi prevalence in prey species (Somov et al. 2019). Sockeye showed significant stimulation of immune and inflammatory genes associated with I. hoferi prevalence.

Similar to I. hoferi, the microsporidian parasite Loma sp. (most likely Loma salmonae) was present at significantly higher prevalence and load in all salmon in the GoA compared with coastal areas (Shaw et al. 2000). This parasite can result in respiratory distress, impaired swimming, and reduced growth rates (Shaw et al. 2000). Transmission is initiated by release of spores from ruptured gill xenomas and is completed by the spores infecting the pillar and endothelial cells of the gills of a new host (Shaw et al. 1998). In the GoA, coho showed significant correlation of gene expression profiles with I. hoferi prevalence.

The bacterium Ca. B. cysticola causes epitheliocystis in gill tissue of salmonids and is associated with proliferative gill inflammation (Toenshoff et al. 2012; Mitchell et al. 2013). This bacterium is commonly encountered in Pacific salmon and was only significantly elevated in prevalence in pink salmon in the GoA (Bass et al. 2017; Teffer et al. 2017). Ca. B. cysticola has been correlated with lower relative weight (Bass et al. submitted), inflammation in coastal Chinook (Wang et al. 2018), and reduced migration success in steelhead (Oncorhynchus mykiss) (Twardek et al. 2019).

Viral encephalopathy and VER was significantly elevated in prevalence in coho salmon. This widespread virus of marine fish and invertebrates is transmitted horizontally, vertically, and trophically (Costa and Thompson 2016). While brain tissue was not included in this screen, in contrast to coastal salmon, the neurotropic VER can also be detected in other tissues by qPCR (Costa and Thompson 2016). We hypothesize that the detection of VER in non-neuronal tissue could reflect a systemic viremia state of recent trophic acquisition, but we likely underestimate both the prevalence and load of this virus in GoA. The relatively high prevalence in coho salmon might reflect their higher trophic level compared with other salmon species encountered in the GoA (Somov et al. 2019).

ENV often causes epizootics in Pacific herring but has recently established as a common coastal virus infecting salmon (Pagowski et al. 2019). It was found in lower prevalence in the GoA, potentially due to the more coastal distribution of Pacific herring not commonly found in the open ocean limiting transmission potential.

The meso/mycetozoea protist Sphaerothecum destruens, transmitted in fresh water in a broad host range of fish, was found at significantly elevated prevalence in sockeye salmon (Gozlan et al. 2009). Infection results in splenomegaly and nephromegaly and causes anemia (Elston et al. 1986). The elevated prevalence in sockeye in the GoA might be a stock of origin, as the GoA has a high proportion of Alaskan origin fish that harbor this infection, compared with the prevalence of British Columbia origin fish in the coastal database.

Ichthyobodo sp. was detected at high prevalence in pink, chum, and coho salmon. This ectoparasite has been shown to be a major factor influencing chum survival at sea in the western Pacific (Urawa 1993; Mizuno et al. 2017).

Interpreting the prevalence data in the GoA compared with coastal British Columbia needs to consider the differences in life stage and season. GoA fish were captured in the middle of their life, whereas coastal salmon from British Columbia were primarily out-migrating post-smolts. As salmon change their diet throughout their life, e.g., increase piscivory, this might impact exposure to trophically transmitted pathogens. The heterogeneity between data sets is especially pronounced for sockeye, pink, and chum salmon that spend most of their life in the open ocean and are only rarely encountered in coastal waters. Coho salmon, on the other hand, are present in coastal waters at all life stages and offer a more robust comparison. Thus, seasonal patterns and fish size or age might influence differences in infectious agent prevalence, specifically for sockeye, pink, and chum.

Myxozoans, commonly observed in coastal environments, have a life cycle that alternates between fish and invertebrate hosts. Most myxozoans, specifically all Parvicapsula spp., showed reduced prevalence in the GoA as invertebrate hosts such as annelids may be limiting (Yokoyama, Grabner and Shirakashi 2012; Somov et al. 2019). As Parvicapsula spp. can reduce visual acuity and have been correlated with increased predation, infected individuals might also be lost from the population (Miller et al. 2014; Nylund et al. 2018). P. kabatai and P. minibicornis both showed stock-specific trends in coho. Tetracapsuloides bryosalmonae, the causative agent of the lethal proliferative kidney disease was absent in the GoA, with infected individuals presumably removed from the population (Sterud et al. 2007).

PSPV, a DNA virus reported in sockeye salmon with unknown pathogenicity, was the highest prevalence virus in the GoA (Miller et al. 2011, 2017; Nekouei et al. 2018). Several novel viruses, Salmovirus and Rhabdo virus, were detected in GoA coho correlating with hypoxia stress and VDD gene expression, as well as a novel Picornavirus in chum (Mordecai et al. 2019, 2020).

The microsporidian P. theridion (syn. Desmozoon lepeophtherii), infects gill tissue but also the sea louse Lepeophtheirus salmonis that may act as a vector (Nylund et al. 2010; Sveen et al. 2012). P. theridion is highly prevalent in coastal salmon in spring and summer but decreases over winter and was only observed in five coho in the GoA (Tucker et al. 2018; Laurin et al. 2019; Bateman et al. 2020).

The bacterium Ca. S. salmonis (Sch), which causes gill impairment, was lower in prevalence in in the GoA (Nylund et al. 2015), as was T. maritimum, the causative agent of mouth rot, presumptively related to the poor outcome of these diseases (Avendaño-Herrera et al. 2006).

In chum and coho salmon, infectious agent and gene expression profiles significantly correlated with size, suggesting that many infectious agents are either shed during maturation or that infected individuals are lost from the population due to mortality. Alternatively, differing prey composition (Losee et al. 2014) or age-dependent mixture such as in chum where Asian origin fish are absent from the first-year age class might explain these trends. Coho salmon showed stock-specific differences in infectious agent profile, with stocks from the contiguous United States showing distinct infectious agent profiles compared to stocks from Northern British Columbia and Alaska.

Changes of the physical environment experienced by salmon at sea based on daily travel rates are negligible (0%–1% on average: Supplementary Material, Table 4, Supplementary Material, Fig. 11) in relation to the speed of gene expression changes that can occur in response to stress in salmonids (Ogura and Ishida 1992; Ogura and Ishida 1995; Houde, Akbarzadeh, et al. 2019). While the abundance of prey items was more spatially variable (3%–28% changes per day on average: Supplementary Material, Table 4, Supplementary Material, Fig. 11), the movement of salmon at sea is not random and salmon are expected to remain in prey rich areas, once found, thus the correlation of gene expression with prey group presence might be stronger than apparent from prey distribution.

Fit-Chip analysis in all species showed large-scale changes in gene expression, specifically from biomarker panels involved in immune response (immune stimulation, inflammation, VDD). In pink, and to a lesser degree in coho salmon primarily caught in warmer waters on the southern border of the survey area, reduced gene expression correlated with warmer temperatures and reduced prey availability. This could be an indicator of higher metabolic demands in malnourished individuals. In chum and sockeye, gene expression correlated positively with temperature, whereas low prey availability still showed a negative correlation. At lower temperatures, gene expression may simply reflect the correlation of metabolic activity with temperature in ectothermic animals. Alternatively, individuals at higher latitudes (i.e., colder waters) were experiencing extremely high abundance of the northern sea nettle Chrysaora melanaster, a large jellyfish (Pakhomov et al. 2019). Thus, temperature might act as a proxy for the impact these large jellyfish had on zooplankton communities thereby affecting lower trophic level salmon in the north of the GoA. Indeed, chum followed by sockeye had the lowest stomach fullness indices (Somov et al. 2019). Individuals with reduced expression of most immune response genes are effectively immunosuppressed. Immunosuppression was correlated inversely with the biomass of the primary prey groups as determined by stomach content in all species except sockeye where small zooplankton had a larger effect than euphausiids, the dominant stomach content of sockeye (Somov et al. 2019). Immunosuppression was strongly correlated with RIB in chum and pink salmon, and to a lesser degree in coho and sockeye. Pink salmon also showed a protective effect of high condition factor that countered immunosuppression and RIB.

Multiple ecological relationships could explain the observed link between energetics (prey availability), immunosuppression, RIB, and temperature. Low prey availability could drive salmon into energetic deficit, to which they respond by suppressing the immune system, a common response to malnutrition in many vertebrates (Latshaw 1991; Lord et al. 1998). Similar observations have been made in steelhead/rainbow trout, where fish exhibit distinct immunity and energetic programs in response to smoltification and migration (Sutherland et al. 2014), as well as in Atlantic salmon where starvation negatively impacted immune response to bacterial infection (Martin et al. 2010). Strikingly, immunosuppression has recently been associated with mortality in Atlantic salmon (Krasnov et al. 2020). Immunosuppression would make salmon more susceptible to pathogens, explaining the elevated infectious agent loads. Immunosuppression could also explain the absence of immune response to pathogens such as Ca. B. cysticola and S. destruens, suggesting that these are opportunistic pathogens with elevated prevalence in immunosuppressed individuals. Since condition factor was inversely correlated with immunosuppression and RIB, “good performance” could have acted protectively, as such individuals are less likely to suffer from energy deficit, thus are immunocompetent and able to fend off infections.

This interpretation is corroborated by field observations, where prey groups showed heterogeneous distributions with little overlap and sockeye and chum salmon exhibited poor feeding condition (Pakhomov et al. 2019; Somov et al. 2019). Specifically in chum, extremely low condition factors individuals were caught where the average water temperature were more than half a degree warmer than their preferred range (Fukuwaka et al. 2007).

Unusually warm temperatures and stratification during the weak 2018–2019 El Niño event—conditions previously hypothesized to disrupt open-ocean food webs and reduce prey availability (Rand 2002; NOAA 2021)—could have driven the observed energy deficits of many salmon in the study area by reducing primary production or altering zooplankton communities. Accordingly, salmon in the survey area were observed to orientate towards structural elements of the water column as well as mixed layer depth, presumably to improve their energetic balance at more favorable environments (Pakhomov et al. 2019; Radchenko et al. 2019). Alternatively, pathogen exposure associated with certain temperature regimes could result in impaired foraging and thereby cause energetic deficits and immunosuppression.

The Fit-Chip technology was developed and validated on the premise of recognizing specific responses based on consistent patterns of coactivation of as few as 7 curated biomarkers (Miller et al. 2017; Houde et al. 2019; Akbarzadeh et al. 2020). However, in the GoA only a subset of any given biomarker panel was co-activated in the first four principal components of gene expression. The observed trends in gene expression were primarily large-scale changes in global gene expression, such as is typical to immunosuppression, rather than responses to specific stressors. One caveat is that this study did not employ known health status controls for different stressors to classify stressor status in individual fish, as these were not available across all four species at the time. We can thus only identify relative differences, rather than classify individuals into specific stressor categories. Refinements of Fit-Chip technology including species-specific stress standards and classification systems are underway.

We presented the first comprehensive overview of the health and condition of Pacific salmon at the end of the winter in the open Eastern Pacific Ocean. We highlight overall trends in pathogen profiles and identify key pathogens present in the open ocean. Further, we find that all species are influenced by energetic constraints correlated with reduced prey availability that was associated with immunosuppression and increased pathogen burden. All species investigated exhibit signs of cumulative effects of stressors, with ocean conditions and prey availability being the primary associated factors. This highlights the impacts a warming ocean could have on winter survival at sea in the face of climate change, specifically in the northern part of the GoA that experienced a large sea surface temperature abnormality in 2019 (Hinch et al. 1995; Miller et al. 2014). Warming, with its downstream effects on salmon energetics, could be especially disruptive in the GoA, where overwintering salmon from both sides of the Pacific basin congregate due to its homogeneous environment (Rand 2002; Beacham et al. 2009; Litzow et al. 2018).

With many wild Pacific salmon populations declining in abundance and productivity, interest in resolving factors that limit salmon survival at sea is strong. Most of what we understand about salmon comes from studies along the coastal margin. The present study provides the first detailed insight into the health and condition of Pacific salmon in the open ocean during the winter. This work will serve as a baseline for future evaluation of the ability of the Northeast Pacific to support salmon populations of North America and Asia.