Network analysis reveals that acute stress exacerbates gene regulatory responses of the gill to seawater in Atlantic salmon

As part of their anadromous life cycle, Atlantic salmon migrate downstream and transition from freshwater (FW) to seawater (SW). To prepare for downstream migration and SW entry, juvenile salmon undergo a developmental period called the parr–smolt transformation. This transformation is triggered by environmental cues including photoperiod and water temperature and consists of several morphological, behavioral and physiological changes (McCormick, 2013). Examples of these changes include changes to the shape and color of the body and fins, increased growth and metabolic rate, loss of territoriality and the acquisition of SW tolerance (McCormick, 2013). The development of SW tolerance involves the ability to withstand direct transfer to SW with minimal osmotic/ionic disturbance and is critical to the performance and survival of salmon smolts once they reach the marine environment.

Upon SW entry, salmon smolts maintain osmotic/ionic homeostasis by employing well-established physiological mechanisms similar to other marine teleosts. In SW, marine teleosts drink, absorb water via the intestine, secrete excess salt via the gill, and produce an isosmotic urine (Grosell, 2010; Hwang et al., 2011; McDonald, 2007). The development of these mechanisms has been extensively studied in salmon smolts and results in part from shifts in the expression and localization of ion transport and water channel proteins in osmoregulatory organs such as the gill, intestine and kidney (McCormick, 2013). For example, during the parr–smolt transformation, protein abundance of the Na⁺/K⁺-ATPase (NKA) (McCormick et al., 2013), the Na-K-Cl cotransporter (NKCC) (Nilsen et al., 2007; Pelis et al., 2001) and the cystic fibrosis transmembrane regulator (CFTR) (Nilsen et al., 2007; Singer et al., 2002) increase in the gill epithelium. Ultimately, these ion transport proteins work together to promote active salt extrusion from the gill and are critical to the SW tolerance of smolts.

The majority of studies examining the physiological mechanisms underlying SW tolerance in Atlantic salmon have focused on the preparatory changes that occur during the parr–smolt transformation while fish are still in FW. However, there are also molecular events that take place in the osmoregulatory organs of smolts after SW exposure that underlie the early phase of SW acclimation. These include further shifts in the gene expression of ion transport proteins (Singer et al., 2002; Sundh et al., 2014; Christensen et al., 2018), as well as the regulation of water channel (aquaporins; Tipsmark et al., 2010a) and tight junction (claudins; Tipsmark et al., 2008, 2010b) proteins within the first 24 h of SW exposure. Together, these studies demonstrate that the remodeling of osmoregulatory tissues continues upon exposure of smolts to SW, playing an important role in the completion of SW acclimation.

The transition from FW to SW represents a physiological challenge for Atlantic salmon smolts, and as such, this life stage is highly sensitive to environmental and anthropogenic stressors (Carey and McCormick, 1998; McCormick et al., 2009; Monette and McCormick, 2008). More specifically, studies have demonstrated a direct link between acute stressors such as increasing water temperature (Vargas-Chacoff et al., 2018), chemical contaminants (McCormick et al., 2005; Nieves-Puigdoller et al., 2007; Waring and Moore, 2004), water chemistry (Monette and McCormick, 2008; Monette et al., 2008, 2010; Nilsen et al., 2010), aquaculture-related handling and transport (Iversen et al., 1998), and the loss of SW tolerance in smolts during the early phase of SW acclimation. From these studies, we know that loss of SW tolerance occurs from impaired endocrine signaling, as well as direct effects on gill ion transport proteins including reductions in NKA activity and NKCC protein abundance, and shifts in the expression patterns of paralogs or ‘isoforms’ of genes that encode for the catalytic α-subunit of NKA. Together, these studies show that acute stressors experienced by smolts prior to or during downstream migration can impair SW tolerance by interfering with well-established hypo-osmoregulatory mechanisms; however, much remains unknown regarding these impacts at the whole-organ level.

Several studies have employed tissue-level transcriptomics to better understand the molecular mechanisms underlying the parr–smolt transformation and the development of SW tolerance in salmonids. Using microarray analyses, Robertson and McCormick (2012) determined that Atlantic salmon smolts exhibit upregulation of genes involved in transcription, protein biosynthesis and folding, electron transport, oxygen transport and sensory perception as compared with the parr life stage. Other studies have determined that the parr–smolt transformation is associated with the regulation of genes related to immunity (Johansson et al., 2016; Sutherland et al., 2014). In Pacific salmon, recent studies have used transcriptomics to identify biomarkers for assessing smolt development (Houde et al., 2019a) and transcriptional responses associated with specific environmental stressors (Houde et al., 2019b). Together, these studies have enhanced our understanding of the parr–smolt transformation and the development of SW tolerance in salmon. However, few studies have utilized genomic tools to examine the regulation of gene networks or pathways underlying these physiological processes. In particular, much remains unknown regarding the impact of acute stressors on the regulation of gene networks involved in SW acclimation.

Our objective was to examine the effects of an acute, aquaculture-related stressor on the dynamics of gene expression in the gills of Atlantic salmon smolts during the early phase of SW acclimation. To examine this, we subjected smolts to handling/confinement followed by a 24 h SW tolerance test. We then combined classic measures of whole-animal physiology with organ-level transcriptomics and gene network analysis to test the hypotheses that: (1) handling/confinement impairs the SW tolerance of smolts, and (2) impaired SW tolerance is associated with shifts in the gene expression response of the gill to SW. Ultimately, the present study sheds light on how stressors encountered prior to or during downstream migration may be impacting the performance of salmon during the early phase of SW acclimation.

Atlantic salmon smolts were obtained from the USDA National Cold Water Marine Aquaculture Center (Franklin, ME, USA) and transferred to the MDI Biological Laboratory (Bar Harbor, ME, USA) in May 2017. Smolts used in this study came from a stock tank containing 101 different families. Prior to the start of the experiment, smolts were held in circular, fiberglass tanks receiving recirculating FW. Tanks were filtered with a biofilter (ProClear Aquatic Systems, Jacksonville, FL, USA), aerated with air stones and maintained at a target temperature of 11°C using a re-circulating chiller system. Fish were exposed to a 12 h:12 h light:dark cycle and fed once a day with 4 mm trout pellets (BioOregon, Westbrook, ME, USA). Fish were held under these conditions for 2 weeks prior to experimentation to allow recovery from transport. Food was withheld for 24 h prior to the start of the experiment. All fish care and protocols were reviewed and approved by the MDI Biological Laboratory Institutional Animal Care and Use Committee (protocol no. c16-03).

To assess the effects of acute stress on SW tolerance, Atlantic salmon smolts (fork length >16 cm) were challenged with 24 h of SW both with and without prior exposure to handling and confinement stress. Just prior to the initiation of the experiment, 10 smolts were sampled directly from FW rearing tanks to serve as the freshwater control group (FW control). At the same time, a second group of 10 smolts was transferred to another circular, fiberglass tank containing artificial SW at a salinity of 35‰, achieved using Instant Ocean Aquarium Salt (Instant Ocean, Spectrum Brands, Blacksburg, VA, USA), and sampled after 24 h to serve as the seawater control group (SW control). A 24 h SW tolerance test as described by Clark (1982) was employed to assess the hypo-osmoregulatory ability of smolts, and as a sampling time-point during the early phase of SW acclimation. A third group of smolts (n=15) was subjected to an acute handling and confinement stress. This stressor was chosen as it has been previously shown to be sub-lethal and to elicit a general stress response in Atlantic salmon smolts (Carey and McCormick, 1998). Fish were captured and held in a net out of water for 20 s and then crowded at an approximate density of 60 kg m⁻³ in a plastic tote (150 l) for 4 h. The plastic tote was covered with a mesh cover and received continuous aeration. During the 4 h confinement period, temperature, dissolved oxygen levels and ammonia levels were continuously monitored. At the end of the 4 h confinement period, seven smolts were sampled directly from the plastic tote to serve as the freshwater stress group (FW stress), and eight smolts were transferred to the 35‰ SW tank and sampled after 24 h to serve as the seawater stress group (SW stress).

All smolts were anesthetized in MS-222 (250 mg l⁻¹ MS-222, Western Chemical, Ferndale, WA, USA), length-measured and weighed. Blood samples were drawn from the caudal vein using a 1 ml ammonium-heparinized syringe within 10 min of initial tank disturbance. Blood samples were centrifuged at 8000 g for 10 min and plasma was aliquoted and stored at −80°C. Gill tissue from the first gill arch was removed above the septum and placed into a 1.5 ml tube containing RNAlater Stabilization Solution (Invitrogen, Carlsbad, CA, USA). Gill samples were stored at 4°C for 24 h and then transferred to −80°C.

To assess the impact of acute stress on whole-animal physiology, we examined plasma osmolality and chloride levels as indicators of osmotic/ionic homeostasis. Plasma glucose was examined in order to confirm that a general stress response was elicited by our handling/confinement stress regime. Plasma osmolality was measured with a vapor pressure osmometer (Wescor 5520 Vapro osmometer, Logan, UT, USA). Plasma chloride was measured by silver titration using a digital chloridometer (Labconco Corp., Kansas City, MO, USA). Plasma glucose was measured by enzymatic coupling with hexokinase and glucose-6-phosphate dehydrogenase (Sigma Diagnostics, 16-UV) as previously described (Slein, 1965).

We employed a fully factorial design of two treatment conditions (control and stress) and two salinities (FW and SW). Two-way ANOVA was used to test for effects of stress, salinity and their interaction using SigmaPlot (version 12.5). Prior to ANOVA, Shapiro–Wilk and Levene's tests were used to determine whether our data met the assumptions of normality and homogeneity of variance, respectively. When these assumptions were not met, data were rank transformed prior to ANOVA. When significant ANOVA effects were observed, a Tukey's post hoc test was used for multiple comparisons.

Total RNA was extracted from gill tissue using the RNAeasy Plus Universal Mini Kit (Qiagen, Germantown, MD, USA) according to the manufacturer's instructions. Clean-up of total RNA was subsequently performed using DNase I (Invitrogen) to eliminate residual genomic DNA. RNA quality of a subset of gill samples (10 of 35) was assessed using a Bioanalyzer (Agilent 2100) on an RNA Nano chip with the total eukaryotic RNA program. RNA integrity numbers ranged from 8.4 to 8.8 for all samples. Purified RNA was sent to the Genomic Sequencing and Analysis Facility at the University of Texas, Austin, for library preparation and sequencing. For this, we employed a 3′ tag-based RNA sequencing approach known as Tag-Seq, which is an efficient alternative to traditional whole mRNA sequencing that has been tested in fish (Lohman et al., 2016) and other vertebrates (Velotta et al., 2020). Libraries were generated from high-quality RNA and each fish (n=35) was individually barcoded with a unique adaptor. All 35 indexed libraries were sequenced on an Illumina HiSeq 2500. Sequencing was conducted to produce 100 base pair reads generating 3–5 million raw pairs of reads per individual. Raw sequences were deposited in the NCBI Short Read Archive (SRP364565). Raw reads from each of the 35 samples were then filtered and trimmed for quality using the FASTX-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Following filtering and trimming, remaining sequencing reads were individually mapped to the Atlantic salmon genome (GCF_000233375.1) using the Burrows–Wheeler alignment tool (BWA) (Li and Durbin, 2010), following the methods of Lohman et al. (2016). Non-primary reads identified as ‘supplementary’ were removed (approximately 4% of mapped reads on average) to avoid counting duplicates or ambiguous mappings. A table of transcript abundances was generated for each gene in the Atlantic salmon genome using featureCounts (Liao et al., 2014). Because genes with low read counts are subject to measurement error (Robinson and Smyth, 2007), we excluded those with less than an average of 10 reads per individual.

Differential gene expression (DE) analysis was performed using the package edgeR (Robinson et al., 2010) in R v 3.6.1 (https://www.r-project.org/). The function calcNormFactors first was used to normalize read counts among all libraries, after which model dispersion for each transcript was estimated separately using the function estimateGLMTagwiseDisp (McCarthy et al., 2012). We tested for differences in transcript abundances between stress, salinity and their interaction using a generalized linear modelling (GLM) approach. GLM log-likelihood ratio tests (LRTs) were used to assess significance after controlling for multiple testing using a genome-wide false discovery rate (FDR) of 0.05 (Benjamini and Hochberg, 1995).

We performed weighted gene correlation network analysis (WGCNA) (Langfelder and Horvath, 2008) to identify the regulation of gene expression associated with responses to stress, salinity and their interaction. WGCNA is a network-based statistical approach that identifies clusters of genes with highly correlated expression profiles (hereafter referred to as modules). We performed this analysis by clustering gene expression profiles of all individuals across all treatment groups; in this way, modules represent clusters of genes whose expression co-varies in response to treatment conditions. This approach has been successfully implemented in studies relating gene expression with whole-organism performance-level phenotypes and physiology (e.g. Velotta et al., 2017, 2020; Campbell-Staton et al., 2020). Prior to performing WGCNA, we normalized raw read counts by library size and log-transformed them using the functions calcNormFactors and cpm, respectively, in edgeR. Module detection was performed using the blockwiseModules function in WGCNA with networkType set to ‘signed’ and otherwise default parameters. Briefly, Pearson correlations of transcript abundance data were calculated between pairs of genes, after which an adjacency matrix was computed by raising the correlation matrix to a soft thresholding power of β=8. Soft thresholding power β is chosen to achieve an approximately scale-free topology, an approach that favors strong correlations (Zhang and Horvath, 2005). We chose a β=8 because it represents the value for which improvement of scale-free topology model fit begins to decrease with increasing thresholding power. A topological overlap measure was computed from the resulting adjacency matrix for each gene pair. Topologically based dissimilarity was then calculated and used as input for average linkage hierarchical clustering in the creation of cluster dendrograms. Modules were identified as branches of the resulting cluster tree using the dynamic tree-cutting method with default parameters (Langfelder and Horvath, 2008). To visualize the co-expression network and quantify the strength of correlations among them, we built an interaction network using the igraph package in R (Csardi and Nepusz, 2006). Connections among modules were calculated as the pairwise Pearson correlation coefficients between module eigengenes. Only significant (P<0.05) connections after FDR correction were retained.

Once modules were defined, we used a multi-step process to associate expression among individuals with variation in measures of whole-animal physiology (plasma osmolality, chloride and glucose). First, we summarized module expression using principal components analysis (PCA) of gene expression profiles (blockwiseModules in WGCNA); because genes within modules are highly correlated by definition, the first principal component (referred to hereafter as the module eigengene) was used to represent module expression (Langfelder and Horvath, 2008). For modules in our data, PC1 explained 27–48% of the variance in gene expression (Table S1). We used module eigengene values to test for associations between module expression and each measure of whole-animal physiology for each module (Pearson correlation; cor function in WGCNA). P-values for the correlation were determined by a Student's asymptotic test (corPvalueStudent in WGCNA). We then conducted ANOVA on rank-transformed module eigengene values in order to test for the effects of stress, salinity and their interaction on module expression. P-values from association tests and ANOVAs were corrected for multiple testing using the FDR method. Finally, we calculated intramodular connectivity values (kME) in WGCNA as the Pearson correlation between an individual gene's expression profile and the module eigengene (Langfelder and Horvath, 2008). kME values are used as a metric of how well individual genes are connected to the modules for which they are assigned.

We performed functional enrichment analysis on all modules using the R package gProfilerR (Reimand et al., 2016). We corrected for multiple testing using gProfilerR's native g:SCS algorithm. We used the list of genes from the filtered gill transcriptomes as a custom background list in all enrichment analyses. We searched for enrichment among terms from the Gene Ontology (GO; The Gene Ontology Consortium 2019) database.

In this study, plasma glucose was used as an indicator of a general stress response in smolts. Two-way ANOVA revealed significant stress (P<0.001) and interaction (P=0.003) effects, but no effect of SW alone (P=0.071) on plasma glucose (Fig. 1A). In FW, exposure to handling and confinement stress caused a 57% increase in plasma glucose in stressed smolts relative to controls, whereas after 24 h in SW, plasma glucose was elevated by 28% in stressed smolts relative to controls (Fig. 1A). A small increase in plasma glucose (19%) was observed with transfer of control smolts to SW (Fig. 1A). Together, these data confirm that our handling and confinement stressor initiated a general stress response, consistent with previous studies in Atlantic salmon smolts (Carey and McCormick, 1998; Iversen et al., 1998, 2005). In addition, sustained elevation of plasma glucose in stressed fish relative to controls after 24 h in SW indicates the persistence of a general stress response during the early phase of SW acclimation despite cessation of the stressor.

Plasma osmolality and chloride were used as indicators of osmotic/ionic homeostasis. Two-way ANOVA revealed a significant effect of stress (P=0.009), SW (P<0.001) and their interaction (P=0.042) on plasma osmolality (Fig. 1B). Plasma osmolality was not significantly different between control and stressed smolts in FW, but was elevated (12%) in stressed smolts relative to controls after 24 h in SW (Fig. 1B). For plasma chloride, two-way ANOVA revealed a significant effect of SW (P<0.001), but not of stress (P=0.163) or their interaction (P=0.100) (Fig. 1C). Similar to plasma osmolality, plasma chloride was not significantly different between control and stressed smolts in FW, but was elevated (12%) in stressed smolts relative to controls after 24 h in SW (Fig. 1C). In control smolts, transfer to SW resulted in smaller increases in plasma osmolality (4.4%; Fig. 1B) and chloride (9.6%; Fig. 1C), consistent with previous observations that fully developed Atlantic salmon smolts exhibit minimal osmotic/ionic disturbance upon transfer to SW (Handeland et al., 1998, 2003). Together, these data indicate that smolts exposed to handling/confinement stress exhibited reductions in the ability to maintain ion/osmoregulatory homeostasis during the first 24 h in SW as compared with controls.

After trimming raw sequence reads, we obtained an average of 2.3 million reads per sample. FastQC results of trimmed sequences showed high sequence quality (Phred score>35), without adaptor contamination or overrepresented sequences (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The proportion of reads mapping to the reference genome was high and generally consistent across samples (average: 91.2±0.33%). The proportion of those counted by featureCounts was necessarily lower (59.6±0.01%), as featureCounts does not count ambiguously mapped reads or reads with low mapping quality (Phred score>20). After filtering out genes that had <10 reads across individuals, 17,851 genes remained.

Of the 17,851 total transcripts identified, 2722 were differentially expressed (1126 downregulated and 1596 upregulated) in stressed smolts as compared with controls (Fig. 2A,B). Exposure to SW alone altered gene expression in the gill, but to a lesser extent than handling and confinement stress, with 1019 genes being differentially expressed (454 downregulated and 565 upregulated) after 24 h in SW as compared with FW smolts (Fig. 2A,B). When handling and confinement was combined with SW, the effect was intermediate, with 1396 genes being differentially expressed (863 downregulated and 553 upregulated) (Fig. 2A,B). Principal components analysis revealed that FW- and SW-exposed fish differentiated strongly along PC1 (which explained 17.4% of variation) whereas control and stress-exposed fish differentiated along PC2 (which explained 13.6% of variation) (Fig. 2C). Differentiation along PC1 between FW and SW was greatest for stressed fish as compared with controls, suggesting that exposure to acute stress exacerbates the differences between FW and SW (Fig. 2C). Together, these data indicate that smolts exposed to acute handling/confinement stress experience large and rapid (within 4 h) shifts in the transcriptional activity of the gill. When both stress and SW were combined, a distinct profile of transcriptional activity emerged, indicating an interactive effect on the gill transcriptome. Full model output for all significantly differentially expressed genes is available in Table S2.

WGCNA was used to identify clusters of co-expressed genes within the gills of smolts that co-vary in response to stress, SW and their interaction. Of the 17,851 gill mRNA transcripts expressed (after filtering), 14,508 genes were successfully assigned into 14 modules designated by number (M1–M14), with the number of genes in each module ranging from 92 to 3127 (Table 1). A full list of gene assignments is available in Table S3. Of the 14 modules identified, 10 (M1, M2, M3, M5, M6, M8, M10, M11, M13 and M14) exhibited both a significant two-way ANOVA effect and a significant correlation to at least one measure of whole-animal physiology (Table 1). Functional enrichment analysis was performed to characterize the putative function(s) of genes in these 10 modules. Of these 10 modules, five (M1, M5, M8, M13 and M14) were significantly enriched for at least one GO Biological Process term (Table S4). To examine correlations between modules, we constructed an interaction network of the 14 modules identified by WGCNA analysis (Fig. 3). Modules that were highly correlated (P<0.001) to each other included M5 and M8, M2 and M11, M13 and M9, and M1, M3, M10 and M14 (Fig. 3). For a full list of module correlations, see Table S5. Together, these modules indicate sets of highly correlated genes with divergent responses to stress, SW and their interaction.

Five modules (M2, M5, M6, M8 and M11) exhibited a significant effect of stress and were correlated to plasma glucose (Table 1). Of these five modules, M2 and M11 were consistent with rapid effects (within 4 h) of handling/confinement stress that were maintained after cessation of the stressor and 24 h in SW (Fig. 4C,D). For M6, handling/confinement stress had a significant impact on gene expression in FW, but this effect was not maintained after 24 h in SW (Fig. 4E). M5 and M8 were consistent with both rapid and transient impacts of handling and confinement stress (Fig. 4A,B). M5 exhibited a negative correlation with plasma glucose (r=−0.40, P=0.03; Table 1), whereas M8 exhibited a positive correlation (r=0.52, P=0.003; Table 1). However, neither module exhibited a significant correlation with plasma osmolality or chloride (Table 1), indicating the expression of genes in response to general stress that is independent of osmotic/ionic homeostasis. As seen in Fig. 4A, genes in M5 were downregulated in response to stress in FW, but returned to control levels after 24 h in SW. A similar pattern was observed for genes in M8; however, in this module, genes were upregulated in response to stress (Fig. 4B), and a significant interaction of stress and SW was observed (P=0.042; Table 1). Together, these data identify gene networks on which acute stress has a large, rapid, but transient effect. However, the interaction effect observed for M8 may indicate that gene expression did not fully recover to control levels after cessation of the stressor and 24 h in SW. Functional enrichment analysis revealed that M5 was enriched for genes primarily involved in RNA biosynthetic and metabolic pathways as well as the regulation of transcription (Table 2,,Table S4). In contrast, M8 was enriched for genes primarily involved in signaling, cell communication, signal transduction, negative regulation of mitogen activated protein kinase (MAPK) activity, and negative regulation of protein serine/threonine kinase activity (Table 2, Table S4).

Patterns of expression on a single module (M13) was consistent with the impacts of SW alone (Table 1). This module exhibited a negative correlation with both plasma osmolality (r=−0.44, P=0.02) and chloride (r=−0.49, P=0.006), but no correlation with plasma glucose (r=−0.15, P=0.50) (Table 1). In this module, gene expression was downregulated in response to SW under both the control and stress conditions (Fig. 5E). Functional enrichment of this module revealed enrichment for genes involved in biological processes including signal transduction, cell communication, regulation of enzyme activity, protein phosphorylation, cytoskeletal organization and semaphorin–plexin signaling (Table 2).

Four modules (M1, M3, M10 and M14) were consistent with combined impacts of stress and SW, as these modules were significantly affected by both stress and SW (Table 1). M1 and M3 exhibited positive correlations with plasma osmolality and chloride, and a positive correlation or a trend towards a positive correlation with plasma glucose, whereas M10 and M14 exhibited negative correlations with all three measures of whole-animal physiology (Table 1). As seen in Fig. 5A, expression of genes in M1 was not affected by acute stress in FW but was upregulated in response to SW. Interestingly, when both stressors were combined, increased mRNA abundance was observed in stressed smolts relative to controls after 24 h in SW. In M3, gene expression was upregulated in response to both stress and SW and exhibited further upregulation with the combination of both stressors (Fig. 5B). A similar response was observed for M10 and M14; however, in these modules, genes were downregulated in response to stress and SW (Fig. 5C,D). Together, these data identify gene networks in which exposure to acute stress has altered gene expression in the gill during the early phase of SW acclimation. Functional enrichment analysis of M1 revealed enrichment for genes involved in multiple biological processes including RNA processing and metabolism, ribosome biogenesis, protein localization, transport and folding. (Table 2, Table S4). M14 was enriched for genes involved in DNA replication (Table 2, Table S4).

To further investigate the impacts of acute stress, seawater and their interaction on gene expression dynamics in the gills of smolts, we calculated the absolute value of intramodular connectivity (kME) in WGCNA to identify hub genes with the strongest correlations (top 5%) to module expression for the modules in Fig. 5. A full list of kME scores for each gene in each module is included in Table S3. We also explored the module location and kME of candidate genes previously shown to be regulated in the gills of Atlantic salmon during the parr–smolt transformation and in response to SW (Table 3). We determined that genes that code for NKCC 1a, FXYD-9 and a V-type H⁺ ATPase subunit exhibited significant correlations to M1 expression (kME=0.79–0.86) (Table 3). In addition, genes coding for NKA β subunit and insulin-like growth factor binding protein-5a exhibited significant correlations to M3 (kME=0.88) and M10 (kME=0.85) expression, respectively (Table 3). Finally, genes coding for tight junction protein ZO-2-like, semaphorin 6D-like and plexin-B2-like exhibited significant correlations to M13 expression (kME=0.82–0.90) (Table 3).

Understanding responses to interacting environmental and anthropogenic stressors is critical to predicting physiological performance. This may be especially true for the smolt life stage of Atlantic salmon, as mortality is estimated to be high for this life stage (Thorstad et al., 2012). In the present study, we examined the interaction of an acute, aquaculture-related stressor and SW exposure on the gill transcriptome of Atlantic salmon smolts. To do this, we utilized a systems-level approach by combining classic measures of osmoregulation/ionoregulation and general stress with organ-level transcriptomics and gene network analysis. We were able to: (1) show that acute handling/confinement stress impairs SW tolerance of smolts; (2) identify gene networks exhibiting rapid responses to acute stress; and (3) determine the influence of acute stress on putative gene regulatory pathways related to early SW acclimation.

Exposure to handling/confinement stress had no detectable impact on osmotic/ionic homeostasis in FW, but impaired the SW tolerance of smolts as indicated by increased osmotic/ionic disturbance in stressed fish relative to controls after 24 h in SW. This result confirms those of previous studies demonstrating that the ability of smolts to maintain osmotic/ionic homeostasis in SW is highly sensitive to environmental and anthropogenic stressors (Iversen et al., 1998; Monette et al., 2008, 2010; Vargas-Chacoff et al., 2018; Waring and Moore, 2004). Given the acute time-course of the stressor used in this study (4 h), we suggest that loss of SW tolerance has resulted in part from the release of catecholamines, and subsequent changes to blood flow and gill perfusion, which are well-established mechanisms used by fish to increase oxygen supply to tissues during acute stress (Nilsson, 2007). Although we did not test this directly, this conclusion is supported by elevated plasma glucose levels in stressed smolts relative to controls in both FW and SW. In fish, elevated plasma glucose levels have been shown to result from the activation of the hypothalamic–pituitary–adrenal axis and the subsequent release of catecholamines during stress (Wendelaar Bonga, 1997). As loss of SW tolerance was associated with elevated plasma osmolality and chloride levels, we suggest that stressed smolts experienced increased osmotic disturbance, which may be an important signal for the regulation of gene expression in the gill along with other proximate and systemic signals (Fiol and Kültz, 2007).

We determined that handling/confinement stress had a large, rapid (within 4 h) and transient impact on gene expression dynamics in the gill, and this was associated with a general stress response in smolts. For M5, we observed the downregulation of gene networks involved in numerous metabolic and biosynthetic processes as well as gene transcription, indicating that acute stress leads to rapid suppression of these biological processes. For M8, acute handling/confinement stress led to the upregulation of gene networks involved in signal transduction, including regulators of MAPK activity and serine/threonine kinase activity. Previous studies have shown that both p38 MAPK and its upstream activators (MAPK kinase kinase) are regulated in the gills of Atlantic salmon in response to environmental stress and are thought to play a role in transcriptional regulation and inflammatory response (Hansen and Jørgensen, 2007; Hansen et al., 2008). In other euryhaline models, the MAPK pathway is a well-known element of osmosensory signal transduction networks that integrate and amplify signals from molecular osmosensors to activate downstream effectors, such as ion transport proteins, in response to osmotic stress (Fiol and Kültz, 2007). Rapid stimulation (within 4 h) of genes in M8 is consistent with their role in signal transduction and aligns with previous observations that regulators of MAPK activity are transiently (peak levels observed between 2 and 8 h after SW transfer) induced in the gill in response to osmotic stress (Evans and Somero, 2008; Fiol et al., 2006). Interestingly, the effect on gene expression in these modules was transient, indicating a stress-specific effect that is alleviated after removal of the stressor. However, the recovery of gene expression in M8 was not fully complete after cessation of the stressor and 24 h in SW, which is consistent with sustained elevations in plasma glucose observed in stressed fish relative to controls after 24 h in SW. Together, these results indicate that the regulation of signal transduction pathways is an immediate or early response of the smolt gill to handling/confinement stress that can persist into the early phase of SW acclimation and that could potentially interfere with mechanisms of osmosensing and signal transduction.

We identified four modules (M1, M3, M10 and M14) with gene expression patterns that were consistent with the combined effects of handling/confinement stress and SW, as these modules were correlated to indices of both general stress and osmotic/ionic homeostasis. In first examining the effect of SW alone in these modules, we observed upregulation of gene networks in M1 involved in RNA processing and metabolism, protein translation, folding, localization and transport, indicating the stimulation of numerous metabolic and biosynthetic processes. Work in Japanese eels has shown that biosynthetic and metabolic regulatory pathways are highly enriched in the gill when fish are acclimated to SW as opposed to FW (Lai et al., 2015), which supports the previous finding of higher metabolic rates of SW-acclimated fish (Zikos et al., 2014). In M1, we also identified the upregulation of several proteins previously shown to be regulated in the gills of Atlantic salmon during the parr–smolt transformation or in response to SW, including putative FXYD (fxyd9b), V-type proton ATPase (Atp6ap1b) and NKCC (nkcc1a) proteins (Table 3).

Surprisingly, none of the four gene modules exhibiting combined effects of stress and SW (M1, M3, M10 and M14) were highly enriched in the biological process of ion transport. M3 and M10 may have been promising candidates for the regulation of ion transport processes as these modules exhibited strong correlations with plasma osmolality and chloride (Table 1) but were not significantly enriched for gene ontology terms. However, hub gene analysis for these modules identified genes related to smolting and SW acclimation in Atlantic salmon. In M3, exposure to SW resulted in the upregulation of a putative sodium/potassium-transporting ATPase subunit β-233-like (AT233) (Table 3). In Atlantic salmon, gill mRNA levels of the β subunit of NKA increase from February to April concurrent with the parr–smolt transformation, suggesting a role of this NKA isoform in the development of SW tolerance in smolts (Seidelin et al., 2001; Nilsen et al., 2007). In M10, we identified the downregulation of a putative insulin-like growth factor-binding protein (igfbp-5a) in response to SW (Table 3). The growth hormone/insulin-like growth factor (IGF) I axis is known to be involved in SW acclimation in fish (McCormick, 2001; Sakamoto and McCormick, 2006), and ifgbp-5a is one of the many IGF binding proteins expressed in the gills of Atlantic salmon smolts (MacQueen et al., 2013). In the smolt gill, igfbp-5a is downregulated upon SW exposure (Breves et al., 2017), consistent with the pattern of gene expression observed in our study, and similar to what has been observed in other species of salmon (Maryoung et al., 2015). In zebrafish, Igfbp-5a is thought to regulate calcium uptake in the gill (Dai et al., 2014); therefore, the downregulation of expression in SW may reflect the increased calcium concentration of SW, and subsequent downregulation of calcium uptake pathways (Breves et al., 2017).

In M13 and M14, exposure to SW alone resulted in the downregulation of gene networks involved in signal transduction, cell communication and DNA replication, consistent with previous transcriptomic studies reporting the regulation of these processes in the gill during SW acclimation in other euryhaline models (Lam et al., 2014; Lai et al., 2015; Gibbons et al., 2017; Evans and Somero, 2008; Whitehead et al., 2012; Vij et al., 2020). In particular, our observation of downregulation of a putative ZO-like tight junction protein (tjb2b) in M13 (Table 3) is consistent with the cellular remodeling of the fish gill during SW acclimation. In SW-acclimated fish, tight junction proteins (e.g. occludin, claudin, ZO-1) are downregulated, leading to a gill epithelium that is ‘leakier’ than that of FW-acclimated fish to facilitate the paracellular movement of Na⁺ down its electrochemical gradient (Chasiotis et al., 2012). Specifically, the tight junction protein ZO-1 is believed to play a dual role in providing structural support to tight junctions by linking occludin and claudin to the cytoskeleton and in signal transduction related to gene expression, cell proliferation and differentiation (Chasiotis et al., 2012). In the gills of FW-acclimated Atlantic salmon, the tight junction protein ZO-1 co-localizes with NKA, suggesting an association with mitochondria-rich cells (MRCs); however, little colocalization is observed in the gills of SW-acclimated fish, indicating negligible ZO-1 associations with SW-type MRCs (Engelund et al., 2012). Together, these results support the conclusion that cell adhesion genes (including candidate tight junction proteins and regulators; Table 3) are downregulated in response to SW.

Interestingly, for M1, M3, M10 and M14, gene expression in response to SW was exacerbated with prior exposure to handling/confinement stress (Fig. 5). This is consistent with the observation that stress exacerbates the differences between FW and SW as shown in Fig. 2C. Specifically, this response indicates that exposure to acute stress prior to SW entry can increase the intensity of SW pathways in the smolt gill. This suggests that acute stress responses may necessitate overcompensation of the normal regulation of hypo-osmoregulation (Gunderson et al., 2016). This overcompensation may be critical to the survival of stressed smolts during acute transfer to SW, as work in other euryhaline models indicates that early activation of regulatory mechanisms at the gene level are crucial to euryhaline teleost fishes exposed to abrupt changes in environmental salinity (Lin et al., 2004). We also suggest that the combined effects of acute stress and SW on the gene expression dynamics of several putative proteins involved in ion regulation in Atlantic salmon may demonstrate an increased investment in this process to compensate for impaired SW tolerance. Future studies are needed to determine whether our observations represent an effect that is specific to aquaculture-related handling/confinement stress or is applicable to other environmental and anthropogenic stressors that impair SW tolerance.

We determined that M13 was highly enriched in genes involved in semaphorin–plexin signaling, which may represent a novel signaling pathway regulated in the gills of smolts during the early phase of SW acclimation. To the best of our knowledge, there is no information about semaphorin–plexin signaling in fish; however, it is well known for its role in the development of the central nervous system, including neuronal axon guidance in other vertebrates (Ohta et al., 1995). Studies have also demonstrated a role for this signaling pathway outside of the central nervous system (see review by Perälä et al., 2012). Specifically, semaphorin signaling molecules along with their cell-surface receptors (plexins and neuropilins) have been found to be both positive and negative mediators of barrier integrity and to regulate tight junction proteins in epithelial and endothelial cells (Treps et al., 2013). Assuming the downstream effectors of semaphorin–plexin signaling are similar in fish, we suggest that this pathway could represent an interesting candidate for the regulation of tight junction proteins including claudins, which play an important role in controlling paracellular permeability and are regulated in the gills of salmon during SW acclimation (Tipsmark et al., 2008).

Our data support the conclusion that acute stressors such as those associated with aquaculture-related handling and transport cause shifts in the gene expression dynamics of the gill during the early phase of SW acclimation in Atlantic salmon smolts. We identified strong links between measures of whole-animal physiology and patterns of organ-level gene expression, which provide valuable insight into the response of the fish gill to multiple stressors. We show that acute handling/confinement stress impairs the SW tolerance of smolts and leads to increased investment in SW-regulatory pathways, which may be critical to restoring osmotic/ionic homeostasis in SW. We hypothesize that increased investment in these processes may represent an energetic trade-off with other fitness-related functions including activity, growth and immunity (Sokolova, 2013). It will be important for future studies to determine: (1) the time-course of transcriptomic recovery following acute stress and SW, (2) whether transcriptomic impacts translate into protein level effects and are linked to long-term impacts on physiological performance, and (3) the impacts of multiple stressors on cell-specific transcriptomes in the gill epithelium. Finally, we confirm previous findings in Atlantic salmon, and identify a novel signaling pathway with a potential role in the regulation of gill remodeling in response to SW.

We would like to thank Brian Peterson and Davin O'Connell from the National Cold Water Marine Aquaculture Center (Franklin, ME, USA) for providing the fish for this study. We would also like to thank Jane Disney, Chris Smith, Karlee Markovich and Cammie Phalen from the Mount Desert Island (MDI) Biological Laboratory (Salisbury Cove, ME, USA) for their assistance with multiple aspects of this project. In addition, we thank Western Connecticut State University (Danbury, CT, USA) undergraduate students Christina Ballatore, Alexandra Helmes, Shane Hoyt, Jake Lipinsky, Kristha Mazariegos, Shannon McFarland, Daniel Suquilanda and Kavya Suresh for their help with fish sampling and tissue processing.