Assessing the effects of genotype-by-environment interaction on epigenetic, transcriptomic, and phenotypic response in a Pacific salmon

Genotype-by-environment (GxE) interactions are non-parallel reaction norms among individuals with different genotypes in response to different environmental conditions. GxE interactions are an extension of phenotypic plasticity and consequently studying such interactions improves our ability to predict effects of different environments on phenotype as well as the ﬁtness of genetically distinct organisms and their capacity to interact with ecosystems. Growth hormone transgenic coho salmon grow much faster than non-transgenics when raised in tank environments


Introduction
Phenotypic plasticity is a characteristic of species or strains in which phenotype (e.g., morphology, physiology, and behavior) is not fixed between environments.Because phenotypic plasticity may confer a fitness advantage by allowing non-genetic modification of traits, its role in modifying evolutionary trajectories through genetic accommodation or assimilation has been a fascinating question in evolution and remains a matter of ongoing discussion (Yeh and Price 2004;Laland et al. 2014;Schlichting and Wund 2014;Levis and Pfennig 2016;Charlesworth et al. 2017;Oostra et al. 2018).Phenotypic plasticity may provide higher adaptive potential than the accumulation of de novo mutations due to its ability to affect multiple individuals rapidly, respond to environment-specific triggers, and harness existing (yet cryptic) genetic variation (Levis and Pfennig 2016).For instance, field studies in the threespine stickleback (Gasterosteus aculeatus) demonstrated the potential of phenotypic plasticity in promoting adaptive radiation of two fish ecotypes through parallel selection of developmental phenotypes (Wund et al. 2008).
Genotype-by-environment (GxE) interaction is a type of phenotypic plasticity influenced by both genotype and environment and where the phenotypic responses of organisms with one genotype are different from those of another genotype between two or more environments-specifically, in a non-parallel fashion.It is thought that gains in selective breeding programs can be hampered when unexpected GxE interactions are not included in genetic evaluations (Mulder 2016).Similarly, in risk assessment of transgenic salmon, GxE interactions [e.g., physiology (Lõhmus et al. 2010), predation (Sundstro ¨m et al. 2007), growth, survival (Vandersteen et al. 2019), and antipredator behavior (Sundstro ¨m et al. 2016)] might lead to mistaken conclusions if non-transgenic and transgenic salmon have non-parallel reaction norms across environments (Devlin et al. 2015).While commercial use of transgenic fish has historically been of concern, transgenic Atlantic salmon (Salmo salar) are now being produced commercially (Waltz 2017) in highly secure land-based culture facilities.If such measures are not applied to the rearing of transgenic fish, these fish may gain access to natural environments.Identifying ecological effects of such an outcome is complicated by GxE interactions affecting phenotype.
The underlying mechanisms of phenotypic plasticity and GxE interactions are of interest as they may provide insight into how the environment interacts with the genome and allows improved predictions of the effects of different environments on phenotype (Nicotra et al. 2010).Many studies have focused on the epigenetic mechanisms underlying phenotypic plasticity [e.g., drought tolerance (Herman and Sultan 2016), fungus in multiple environments (Kronholm et al. 2016), migration (Baerwald et al. 2016), growth in various sugar concentrations (Herrera et al. 2012), and effect of rearing environment (Luyer et al. 2017)].Potentially other mechanisms influencing gene expression could underlie phenotypic plasticity and GxE interactions (Smith and Kruglyak 2008;Gutierrez et al. 2009).
Epigenetics is a broad field that can refer to multiple and varied pathways: from DNA methylation to self-regulating transcription factors (reviewed in Martin and Fry 2018;Cavalli and Heard 2019).Epigenetic gene regulation differs from general gene regulation by duration, where effects might extend long after the initial signal has faded (see Cavalli and Heard 2019 for other definitions).Plasticity is related to epigenetics in that it is the ability of regulatory mechanisms (e.g., transcription factors) to be influenced by environmental differences to change or reverse epigenetic factors already in place and governing cellular characteristics (Cavalli and Heard 2019).
For environmentally induced epigenetic changes, DNA methylation is the best known and studied epigenetic marker (Martin and Fry 2018).Methylation occurs when a methyl group is added to a cytosine by a DNA methyltransferase that has been triggered by environmental signals (e.g., chemicals interacting with DNA methyltransferases, availability of the methyl donor SAM, or the availability of cytosines due to transcription factor binding or the presence of trimethylated histone H3 Lys4 near specific genes; Martin and Fry 2018;Greenberg and Bourc'his 2019 ).In humans, methylation appears to play a role in regulating genes like insulin-like growth factor 2 (i.e., less DNA methylation after exposure to famine) and other genes depending on environmental exposure (Heijmans et al. 2008;Martin and Fry 2018).DNA methylation has typically been associated with transcriptional repression by decreasing transcription factor binding, increasing heterochromatin formation, and by recruiting proteins that silence gene transcription (Greenberg and Bourc'his 2019).DNA methylation is likely more complex than this simple model.For example, in human phagocytic cells, gene transcription appeared to be independent of DNA methylation (Pacis et al. 2019).Although this may be due to detection limits and the statistics used in identifying differentially methylated regions (including which methylation mark is being monitored) (Pacis et al. 2019), it is also likely gene expression changes are regulated through many pathways.
Transgenic salmon have dramatically faster growth due to the integration of the growth hormone transgene construct at a single genomic locus (Devlin et al. 1994;Phillips and Devlin 2010).In the wild, non-transgenic salmonids possess two growth hormone genes (GH1 and GH2) that influence many biological processes (e.g., growth, osmoregulation, and immune function) and are central to a complex signaling system that integrates information from multiple sources (Bjo ¨rnsson et al. 2002).Central figures in relaying GH signaling in salmonids are the GH-binding proteins (GHBP; regulates circulating GH), growth hormone receptor (GHR; highest density in the liver), and insulin-like growth factor 1 (IGF1; gene expression stimulated by GH and primarily produced in the liver) (Bjo ¨rnsson et al. 2002).The transgenic coho salmon in this study have an extra and constitutively active versions of the growth hormone gene (Devlin et al. 1994) causing growth hormone transgene to be expressed in all tissues and is known to have strong effects on gene expression, growth, and physiology (Rise et al. 2006;Devlin et al. 2009Devlin et al. , 2015)).
In previous studies (Sundstro ¨m et al. 2007(Sundstro ¨m et al. , 2016;;Vandersteen et al. 2019), GxE interactions were observed for various traits when non-transgenic and transgenic coho salmon were reared both in stream and in tank environments.Specifically, the elevated growth of transgenic individuals seen in simple tank environments was not observed in the naturalized and complex stream environment where non-transgenic and transgenic coho salmon grew at the same rate (Sundstro ¨m et al. 2007).The plasticity and GxE interaction effects observed in GH transgenic salmon provide a useful model system to better understand epigenetic mechanisms affecting phenotype in vertebrates.We sought to identify underlying gene expression differences related to GxE interactions and to assess if gene expression GxE interactions were influenced by differential methylation.

Materials and methods
Full details of the materials and methods used in this study can be found in File S1.Briefly, growth hormone transgenic and nontransgenic coho salmon (Sundstro ¨m et al. 2007) were reared in tank or simulated stream environments.Liver tissues were collected from two sets of six individuals from each environment and genotype and each group of six was pooled as a single replicate (two replicates for each environment and genotype).RNAseq and whole-genome bisulfite sequencing was performed on each replicate (NCBI SRA accessions: SRX9236012-SRX9236027).The analyses of the RNA-seq and whole-genome bisulfite sequencing data included pairwise comparisons and a genotypeby-environment analysis.Methylation data was analyzed based on promoter region, gene body region, and intergenic sites (scripts can be found in File S2).Gene ontology analyses were used to identify enriched categories of genes.

Data availability
Sequence data were deposited to the NCBI's sequence read archive available through the BioProject ID: PRJNA667073, and R/ python scripts are available as supplementary files (File S2-GSA figshare).All supplemental files can be accessed through figshare: https://doi.org/10.25387/g3.13542158.

Body weight and specific growth rate
Body weights did not reflect growth rate differences among environments or genotypes as salmon groups were size-matched to avoid confounding effects of developmental stage which occurs when analyzing strains with different growth rates sampled at the same age (Figure 1A).A significant interaction between environment and genotype (P < 2e-16) for specific growth rate (SGR) was identified using a two-way ANOVA test (Figure 1A).The average non-transgenic SGR was increased relative to the average transgenic SGR in the stream environment but decreased in the tank environment (Figure 1B).

RNA-seq and WGBS sequence data
RNA-seq paired-end counts ranged from 36.6 M to 59.2 M per pooled sample (NCBI BioProject: PRJNA667073).After trimming, a mean of 42.9 6 2.5 M reads were retained per pooled sample.Read alignment rates averaged 89 6 3.8% for properly-paired reads mapping to the coho reference genome.WGBS paired-end sequence counts ranged from 481.9 M to 512.5 M per pooled sample.After trimming, a mean of 441.6 6 4.1 M reads were retained per pooled sample.Overall, mapping rates were 45.8 6 0.23% for uniquely and paired reads [around 17Â coverage assuming an estimated 2.4 Gbp genome, which is likely a higher coverage than necessary for analysis (Ziller et al. 2015)].

RNA-seq analysis
Transcription of genes in the liver was greatly influenced by the rearing environment (Figure 2A, Table 1).Replicate groups clustered closely indicating shared responses to both genotype and environment.Only 382 DEGs were identified between transgenic and non-transgenic fish in the stream environment compared to 3134 differentially expressed genes (DEGs) in a tank environment (Figure 2B, Table 1).This difference between environments was consistent with the growth rate analysis, where transgenic and non-transgenic fish had similar SGRs in a stream environment, but more divergent SGRs in a tank environment (Figure 1A).Three enriched biological process GO terms were identified between transgenic and non-transgenic salmon in the stream environment: complement activation, 'complement activation, alternative pathway', and organic hydroxy compound catabolic process.In the tank environment, 304 biological process enriched GO terms were identified (Supplementary File S3, Figure S1), many of which were shared with GxE enriched gene categories (see below).There were 759 shared DEGs between this comparison and the GxE comparison (84% of the GxE DEGs-see below).
Both non-transgenic and transgenic salmon had a similar number of liver DEGs between the stream and tank environments but only about 50% of those DEG were shared between both groups (Figure 2B, Table 1, non-transgenic: 5826, transgenic: 6216, shared: 3083).There were 779 enriched biological process GO terms identified between environments for non-transgenic  salmon and 767 for transgenic salmon, with 452 shared (Supplementary File S3, Figure S2).Major themes from both sets of GO terms include: metabolism, ribosome biogenesis, stress, immune system, cell death, chromatin organization, and gene expression regulation (including epigenetic regulation).
From the ANOVA-style analysis (any gene that was differentially expressed between any of the comparisons), 9615 DEGs were identified and 901 DEGs were classified from the GxE analysis (Figure 3A, Supplementary File S3).There were 311 enriched biological processes, 7 cellular components, and 37 molecular function GO terms identified from the GxE analysis (Figure 3B, Supplementary File S3).The GO term with the lowest P-value was response to hormone and is consistent with the growth hormone transgene physiological alteration observed in SGR (Sundstro ¨m et al. 2007).Many of the other lowest P-value enriched GO terms were responses to various stimuli (e.g., lipid, insulin, stress, and nutrient levels), signaling (e.g., insulin receptor signaling pathway, JAK pathway signal transduction adaptor activity, and interleukin-7-mediated signaling pathway), and regulation (e.g., regulation of primary metabolic process, regulation of complement activation, and regulation of hormone secretion) were other main themes found in these GO terms (File S3).
To better understand GxE transcription interactions, the 901 GxE DEGS were clustered based on transcript expression patterns (Figure 4, A and B).Out of 24 identified clusters, enriched GO terms could be classified for 11 of them (Figure 4C).For clusters 1-3, transgenic transcript expression tended to increase from a stream to tank environment, while expression decreased for nontransgenic fish between environments (Figure 4C).The lowest Pvalue categories for each of these clusters were positive regulation of vitamin D receptor signaling pathway, regulation of cellular metabolic process, and negative regulation of interleukin-2 secretion in increasing order (File S3).
The GO category with the lowest P-value from cluster 8, ribosome biogenesis, is suggestive that transcription of genes involved in ribosome biogenesis increases in both transgenic and nontransgenic salmon in tank environments, but to a greater extent in transgenic fish (Figure 4C, Supplementary File S3).This pattern of transcription mirrors SGR between environments.Transgenic coho salmon have decreased transcription of genes related to clusters 9-10 in a tank environment.GO categories with the lowest P-value from these clusters include, in increasing cluster order, digestion and complement activation (File S3).Transgenic fish displayed an opposite transcription pattern for these categories.Transcription in cluster 11 increases from stream to tank environments in transgenic salmon and decreases in non-transgenic salmon (similar to clusters 1-3).The top enriched GO category for this cluster was regulation of cardiac chamber formation.
Clusters 13, 21, and 22 all have similar gene transcription patterns with increased transcription in a tank environment for nontransgenic fish and decreased for transgenic fish.Top enriched GO terms for these clusters were acute inflammatory response to antigenic stimulus, positive regulation of cell death, and glucose metabolic process, respectively (File S3).Hepatic gene transcription decreased for genes in cluster 20 in a tank environment for transgenic and non-transgenic salmon, but the drop in expression was   greater for non-transgenic salmon.The top enriched GO term for this cluster was response to cesium ions (File S3).

WGBS analysis
Methylation patterns differed between both promoter regions and intergenic sites compared to gene regions (Figure 5, Supplementary Figures S3 and S4).As was the case for RNA-seq analysis, replicate groups clustered closely.Most of the changes in hepatic promoter methylation and intergenic site methylation occurred between transgenic and non-transgenic salmon in both stream and tank environments (Figure 5B).There were 4 out of 382 DEGs (1%) that were also differentially methylated between non-transgenic and transgenic salmon in the stream environment and 34 out of 3134 DEGS (1%) in the tank environment.Methylation at 18 promoter regions displayed GxE interaction (Table 2, Supplementary Figure S5).In total, 823 promoter regions were differentially methylated in the ANOVA-style analysis (differentially methylated in any of the comparisons).No enriched GO terms were identified for any promoter methylation comparison.GxE methylation patterns were identified at 485 intergenic sites, and 76050 intergenic sites were found to be differentially methylated in the ANOVA-style analysis (File S4).
In gene regions, most changes in methylation were observed between environments rather than between genotypes as seen in promoter regions and intergenic sites (Figure 5B).More than 95% of the observed changes were from increased methylation of genes in the stream environment.Out of the 6216 DEGs found between environments for transgenic salmon, 221 (3.6%) were also differentially methylated, while 222 out of the 5826 DEGs (3.8%) for the non-transgenic salmon were differentially methylated.Common enriched GO terms found between stream and tank environments for both transgenic and non-transgenic comparisons include: behavior, locomotion, nervous system development, and homophillic cell adhesion via plasma membrane adhesion molecules (Figure 6, Supplementary File S4).Nine GO terms from these two comparisons were shared with those found in the GxE gene expression analysis, but these terms were from broad categories (e.g., cell differentiation, cellular developmental process, and developmental process).There were 63 gene regions with GxE interactions (File S4).No enriched GxE GO terms were identified, but there were several important transcription factors influencing growth identified by this analysis including: neurogenic differentiation factor 1-like, transcription factor jun-B-like, and cyclic AMPdependent transcription factor ATF-1-like (Supplementary Figure S6).Only the neurogenic differentiation factor 1-like gene had methylation patterns expected to influence gene expression (i.e., an increase in methylation corresponded to a decrease in gene expression between environments and a decrease in methylation corresponded to an increase in gene expression) (Supplementary Figure S6).

Discussion
The goal of this study was to better understand how vertebrate phenotypes can be affected by genotype and environment and their interactions at a genomic level.We have utilized a GH From the RNA-seq data, we identified genes with genotype, environmental, and GxE interaction transcription patterns in the liver.Clustering genes by GxE interaction helped us to identify which biological processes changed in a similar fashion as seen for growth.To better understand how genes with GxE interactions were regulated, whole-genome bisulfite sequencing was performed on the same samples and tissue.Although we observed that many genes had increased methylation in the stream environment compared to the tank environment and reduced expression for those with increased methylation overall, we found little overlap with changes in methylation and changes in gene expression.Each of these observations are described in more detail below.

Specific growth rate
Both non-transgenic and transgenic salmon grew faster in the tank environment than in a stream environment.Environmental, physiological, and behavioral factors likely contribute to the twofold increase in SGR in tanks [e.g., availability of food, nutritional capabilities, and suppression of predator avoidance (Devlin et al. 2015)].To uncover which genes influence or underlie these responses, we searched for genes with liver transcription profiles with complex responses including GxE interactions.

RNA-seq genotype differences
In the stream environment, only three enriched GO terms (complement activation; complement activation, alternative pathway; and organic hydroxy compound catabolic process among 382 DEGs) were identified from the comparison between non-transgenic and transgenic salmon liver tissue.Two of these were related to complement activation, suggesting that some ) for all regions or sites in the liver of transgenic and non-transgenic coho salmon reared in either stream or tank environments (expressed in log fold change, logFC).The first MDS plot was generated from all promoter regions, the second was generated from all gene regions, and the last from all methylation sites outside of gene and promoter regions.(B) Venn diagrams of overlapping promoters, genes, or intergenic sites found between different comparisons.Below the promoter and gene Venn diagrams are lines showing the number of overlapping DM regions between promoter and gene based analyses.

NCBI Accession
differences between liver transcription of transgenic and nontransgenic salmon in the stream environment could be due to environment specific pathogens (Thorgersen et al. 2019).In the tank environment, differences in B cell and T cell regulation were noted between transgenic and non-transgenic salmon, but not with the complement activation GO term.As these immunerelated differences were specific to only one of the environments, they represent GxE interactions and were also detected in the GxE analysis.
Recent studies found that GH transgenic salmon have altered immune responses to pathogen mimics in tank environments (Alzaid et al. 2018;Kim et al. 2019).Taken together with the current study, transgenic, and non-transgenic salmon appear to respond differently to pathogen exposure in multiple environments.In the current study, the transgene appeared to influence immune system differences while the environment appeared to influence which aspect of the immune-system (complement activation vs T cell or B cell regulation) was different between transgenic and non-transgenic salmon.The environmental differences could be driven by different pathogen communities, but it could also be as simple as temperature differences (Ignatz et al. 2020).Alzaid et al. (2018) hypothesized that the immune system may modulate growth as a trade-off.Further research into these immune differences may improve cultivation practices of farmed salmon and risk assessments of transgenic salmon as immune function is important to both objectives.
Many of the enriched GO terms identified in the comparison between transgenic and non-transgenic coho salmon in the tank environment were also shared with the GxE comparison (115 shared out of the 304 in the tank environment).Most of the remaining GO terms, without exact matches, shared similar biological functions (e.g., regulation of mitotic cell cycle vs regulation of mitotic cell cycle phase transition).For this reason, discussion of these GO terms is reserved for the GxE section below.It is important to note here, however, that these differences were not seen in the stream environment.This suggests that controlled studies in tank environments may not accurately represent results relevant to stream or wild environments, and, as such, the best risk assessments would incorporate data from multiple environments.

RNA-seq environmental differences
An unanticipated finding from comparing stream and tank environments was that there might be an influence of food availability in the stream environment, as the cellular response to starvation enriched GO term was identified between environments for both the transgenic and non-transgenic salmon.This result was surprising because transgenic and non-transgenic salmon both had positive growth in stream and tank environments.Also, transgenic and non-transgenic salmon were fed to excess in both environments (adding more natural food to this experimental stream system does not increase growth rates; Vandersteen et al. 2019).We observed that in a stream environment both genotypes had reduced growth relative to the tank environment, but both environments still had positive growth overall.Some questions arise from these observations, e.g., are the fish in the stream environment responding specifically to reduced food supply, and are the cellular mechanisms initiated by lower levels of available food also activated by other conditions or under submaximal feed regimes where fish must compete for available natural prey?
There are environmental and behavioral differences that result in reduced nutrition between stream and tank environments [e.g., the nutrition of the prey in the stream environment may be different than feed pellets causing decreased growth rate (Holm et al. 1984;Vandersteen et al. 2019)].It is possible that submaximal feed levels in the stream environment initiated cellular responses similar to starvation if these are more general purpose cellular responses and do not strictly respond to complete nutrient deprivation.This explanation fits with the observation of reduced growth in the stream environment, the cellular responses to starvation GO term, and the feeding regime used in this study.
Many other shared differences of non-transgenic and transgenic salmon were observed between stream and tank environments.They include enriched GO categories related to circadian rhythm, metabolism, protein folding, autophagy, cellular response to stress, chromatin silencing, demethylation, and Figure 6 Enriched GO terms identified from genes with differential methylation between stream and tank environments.GO terms complexity was reduced using Revigo software (Supek et al. 2011), and GO terms with labels were chosen to maximize the shared GO terms and based on ease of interpretation or P-value.The x and y-axes are coordinates from multidimensional scaling of GO term semantic similarities into two dimensions.(A) Enriched GO terms identified between the stream and tank non-transgenic comparison of gene body methylation.(B) Enriched GO terms identified between the stream and tank transgenic comparison of gene body methylation.regulation of gene expression (epigenetic) among many others.While these GO terms may be informative regarding trait differences between fish reared in tank and enriched stream environments, they do not offer obvious insight into growth (the trait under investigation).

RNA-seq GxE interactions
From the 901 DEGs with GxE transcription interactions in the liver and corresponding enriched GO terms, several inferences can be made regarding the physiological and gene transcription responses regulating the SGR GxE interaction.First, the enriched GO term with the lowest P-value, response to hormone, supports that the growth hormone transgene was a main genetic influence in the SGR GxE interaction.Other enriched GO terms further support the influence of the growth hormone in the SGR GxE interaction including: JAK pathway signal transduction adaptor activity, ERK1 and ERK2 cascade, and insulin receptor signaling pathway (insulin-like growth factor I was one of the 901 GxE DEGs) as expected from the known effects of growth hormone (Dehkhoda et al. 2018).This result was expected since salmon with similar genetic backgrounds were used-with the exception of the transgene.Second, several other enriched GxE interaction enriched GO term categories (e.g., metabolism, glycogen metabolism, gluconeogenesis, response to insulin, response to cAMP, response to glucocorticoid, negative regulation of growth, and response to nutrient levels) implicate genes involved in nutrition as downstream drivers of the SGR GxE interaction.Several of these enriched GO terms were found within clusters 21 and 22, which helps us to better understand gene expression patterns associated with these GO terms.In clusters 21 and 22, gene transcription has an inverse relationship to SGR patterns between environments (i.e., gene transcription is lower in the stream environment and higher in the tank environment for the non-transgenic salmon compared to transgenic salmon).One explanation for these expression patterns and the enriched GO terms is that glycogen metabolism and gluconeogenesis may be occurring to a greater extent in liver tissues of transgenic salmon in a stream environment and in non-transgenic salmon in the tank environment.
In humans, glycogenolysis in the liver provides glucose to the rest of the body during short fasting periods and gluconeogenesis provides glucose to the rest of the body for longer starvation periods or with increased stress (Zhang et al. 2019).In the current study, growth was positive in both stream and tank environments, and with feed in excess in the tank environment, starvation seems unlikely (see above for discussion on starvation).The enriched GO terms cellular response to glucocotricoid stimulus and cellular response to stress in clusters 21 and 22 are indicative of a stress response, which may also promote gluconeogenesis and glycogen storage (Kuo et al. 2015).Stress may influence growth by reducing attempts to feed or by reduced feed absorption (Van Weerd and Komen 1998).Stress has also been shown to cause cortisol and glucose GxE interactions in other fish species (Van Weerd and Komen 1998).Stress responses may help explain the difference in growth data seen between transgenic and nontransgenic salmon from other studies under various conditions and in different environments (Vandersteen et al. 2019).Similarly, cortisol and glucose responses can vary based on stressor, genotype, and early life history (Van Weerd and Komen 1998) making stress a good candidate for the growth discrepancy observed in previous studies.We note that GH transgenic coho salmon do show altered carbohydrate (complex and glucose) metabolism, but GxE interactions have not been assessed in these studies (Higgs et al. 2009;Panserat et al. 2014).
With regards to the GxE enriched GO terms, liver protein synthesis has an association with SGR.The gene transcription patterns in cluster 8 mirror the SGR patterns, and this cluster of genes is enriched for genes involved in ribosome biogenesis.Growth and ribosome biogenesis is intimately linked as protein production is necessary for cell growth (Lempia ¨inen and Shore 2009).
From the enriched GO term and clustering analyses, we were able to link SGR GxE patterns to gene expression patterns of the genes downstream of the growth hormone, stress-related genes, and genes responsible for ribosomal biogenesis.We also identified several other categories of genes with GxE transcription interactions that may not have a direct influence on SGR, but may be important indirectly or for other unmeasured traits.These include genes with immune function (e.g., complement activation, regulation of immune system process, response to interleukin-7, and defense response), tissue regeneration (e.g.animal organ regeneration and liver regeneration), vitamin metabolism (e.g., vitamin D3 metabolic process, vitamin D catabolic process, and response to vitamin B2), and pigment formation (e.g., pigment metabolic process).Future research will be needed to better understand why these gene categories also showed GxE interactions.

WGBS genotype differences
While there were extensive differences between non-transgenic and transgenic methylation patterns in promoter regions and intergenic sites, there were no enriched GO terms identified for the promoter regions.There is likely important biological significance and possibly overlap with the SGR phenotype due to these differences, but the evidence suggests they underlie only a fraction of gene expression differences in the liver between transgenic and non-transgenic salmon.For example, gene transcription is markedly different between non-transgenic and transgenic salmon in the tank environment (3134 DEGs), but similar in a stream environment (382 DEGs).Meanwhile, there are similar levels of differentially methylated promoters and intergenic sites between non-transgenic and transgenic salmon in tank (366 differentially methylated promoters and 31,016 intergenic sites, File S5) and stream environments (326 differentially methylated promoters and 29,194 intergenic sites).While the overall increased methylation and reduced gene expression (see below) in the stream environment for these genes with increased methylation could also be responsible for the few DEGs in the stream environment, we did not observe a high proportion of overlapping DEGs between environments and genes that were differentially methylated (e.g., 221 differentially methylated gene regions shared out of 6216 DEGs between environments for transgenic salmon, File S4).One explanation for the divergence between liver gene transcription and methylation data is that the data reflects unobserved factors at the level of the organ or whole organism not accounted for with single tissue data.This hypothesis is supported by gene region methylation patterns seen between environments (see below).

WGBS environmental differences
Distinctive liver methylation patterns for gene regions were observed between stream and tank-reared salmon for both non-transgenic and transgenic salmon.Many of the enriched biological process GO terms between environments were shared between non-transgenic and transgenic salmon (341 non-transgenic, 215 transgenic, 154 shared).Surprisingly, the biological process GO terms between environments were mostly unrelated to known liver biological processes (e.g., learning or memory, muscle contraction, and thigmotaxis).More than 95% of the gene regions involved in these biological processes had increased methylation in the stream environment.Because theses GO terms are likely important for other organs rather than for liver tissues and the genes from these terms have increased methylation, we suggest that they are under strict gene transcription repression in the liver due to environmental factors [e.g., nutritional differences between prey and hatchery feed (Holm et al. 1984;Vandersteen et al. 2019)].The average counts per million (CPM) value for all genes and samples after filtering was 38.4 (n ¼ 26,084, stdev: 493.4), whereas for the subset of genes found to have significantly increased methylation in a stream environment, the average CPM was 33.8 (n ¼ 1262, stdev: 523.1) for the non-transgenic comparison and 11.4 (n ¼ 1250, stdev: 104.3) for the transgenic comparison.In a previous study, hatchery coho salmon (analogous to tank-reared in this study) had increased methylation of most differentially methylated regions in muscle tissue compared to natural salmon (analogous to stream-reared) (Luyer et al. 2017).While the pattern of methylation was the opposite for this study, it is interesting that the trend for both tissues was unidirectional.

WGBS GxE interactions
There were 18 promoter regions and 63 gene regions identified with GxE methylation patterns.None of these promoter or gene regions overlapped with genes with GxE gene expression patterns, and only three promoter and ten gene regions with methylation GxE patterns were found to be differentially expressed in an ANOVA-style analysis of gene expression data.Of the 901 genes with GxE patterns, 23 promoters and 113 gene regions were found to have differential methylation in an ANOVA-style analysis.Both observations lead us to conclude that the majority of GxE transcription was not influenced directly by methylation and that other mechanisms influencing gene transcription were likely responsible.The influence of long-distance methylation and methylation of hub genes involved in large networks may still play important regulatory roles in the observed GxE transcription, but were not examined in this study.Another possibility not explored in this study, but may add additional insight is the influence of methylation position (e.g., could changing the position of methyl groups within the promoter region change gene expression while not appearing to be significantly different).As mentioned in the Introduction, gene expression and DNA methylation may not have a direct relationship, as exemplified in human phagocytes (Pacis et al. 2019), and subtle or complex relationships may not be accounted for in the type of analysis used.However, a disconnect between methylation and gene expression has also been observed in other studies, where gene expression differences were not associated with methylation differences for the majority of genes (Cunningham et al. 2019;Natri et al. 2020).These observations do not suggest that DNA methylation did not play an important role in gene expression regulation in general, but that the gene expression differences were not obviously linked to changes in methylation.For example, in the promoter region there was two clusters of genes (Supplementary Figure S3), one with a peak with around 25% average methylation and another with 90% average methylation.The first cluster had a peak at around 10 counts per million (CPM, measure of transcription), while the second peak was below 1 CPM.These clusters are consistent with expectations of decreased gene expression with increased promoter methylation and was seen in the majority of genes.
Methylation of DNA coding transcription factors may still play a significant role in the observed SGR GxE interaction.For example, two transcription factors (transcription factor jun-B-like and neurogenic differentiation factor 1-like) with GxE methylation patterns in gene regions may both play roles in growth (Andrali et al. 2008;Raffaello et al. 2010).Transcription factor jun-B plays an essential role in skeletal muscle mass maintenance, with increased expression associated with increased muscle mass and decreased expression with atrophy (Raffaello et al. 2010).Transcription of this factor is congruent with the SGR GxE interaction (Supplementary Figure S6), but did not reach significance in the GxE analysis.Methylation differences do not correlate with the expression pattern with this gene.Neurogenic differentiation factor 1-like regulates insulin gene expression in the pancreas and liver (Andrali et al. 2008).While this gene did not have significant GxE gene expression, it did have a gene expression pattern that was the opposite of the SGR GxE profile.If DNA methylation of gene or promoter regions in the liver underlies the SGR GxE profile, it likely does so through a gene that can regulate many other genes (i.e., a gene hub) like these transcription factors because there is little evidence for methylation differences that could explain the 901 GxE DEGs otherwise.More likely, DNA methylation and other regulators of gene expression play an elaborate role in the observed liver GxE gene expression profile.

Conclusions
Liver gene transcription data in transgenic and non-transgenic salmon pointed to genes related to stress as candidates underlying the GxE pattern of growth in this study.For the model system used here, if transgenic and non-transgenic coho salmon respond to environments differently in terms of physiology or behavior based on stress, incorporating components such as stress proxies into predictive models may, for example, increase accuracy for how a transgenic salmon will respond to environmental conditions and thus improve our understanding of the risks transgenic salmon might pose to wild populations.Interestingly, DNA methylation appeared to be independently regulated between promoter and gene regions in the liver and did not seem to play a direct role in the gene transcription responses observed in this study.While not directly influencing gene expression differences seen in the liver, this does not preclude DNA methylation as a common mechanism for gene regulation.It does suggest, however, that other gene regulatory systems, potentially originating in other organs, are responsible for gene expression differences observed in the liver.The present study has shown the remarkable responses in gene transcription that arise from environmental, genetic, and GxE influences, and has revealed underlying complex genomic mechanisms affecting plasticity and GxE interactions influencing both morphological and physiological phenotypes.

Figure 1
Figure 1 Genotype-by-environment influence on weight and specific growth rate of coho salmon.(A) Average coho salmon size per group were matched to avoid confounding effects of developmental stage.No variables were significantly different for weight.(B)The specific growth rate for each group is shown (mean 6 SEM and individual specific growth rate).Environment, genotype, and the interaction of both were all significant with P < 2e-16, P ¼ 8.59e-15, and P < 2e-16, respectively.

Figure 2
Figure 2 Multidimensional scaling plot (MDS) of gene expression data and a Venn diagram of overlapping differentially expressed genes (DEGs).(A) A MDS plot of transcript counts for all genes expressed in the liver of transgenic and non-transgenic coho salmon reared in either stream or tank environments (as expressed in log fold change, logFC).Each point in the MDS plots represent six pooled individuals.(B) Overlapping DEGs between comparisons.

Figure 3
Figure 3 Genotype-by-environment interactions of liver mRNA transcription in a stream and tank environment.(A) A heatmap of transcript counts per million (CPM) values of 901 DEGs with GxE interactions.The NT symbol represents non-transgenic and T represents transgenic.For each experimental condition there were two replicates of pooled individuals.(B) Enriched GO terms of the genes with GxE interactions after complexity reduction [i.e., cluster representatives are shown (Supek et al. 2011)] of the 311 biological process GO terms originally identified.The x and y-axes are coordinates from multidimensional scaling of GO term semantic similarities into two dimensions.The size and color of each GO term bubble reflects the P-value, with larger and blue bubbles with the lowest P-values.

Figure 4
Figure 4 Genotype-by-environment interactions of liver mRNA-transcription clustering.(A) A heatmap of Spearman correlation coefficients with genes clustered with similar gene transcription.Each of these genes had transcript GxE interaction.The NT symbol represents non-transgenic and the T represents transgenic.(B) Two example clusters are shown of the 24 identified.Colors of the heatmap are based on transcript counts per million (CPM) values.(C) Illustrations of expression patterns for each cluster found to have enriched GO terms (does not reflect actual expression values).Some clusters were grouped in this figure, but each were individually analyzed for GO term enrichment.Boxes above graphs show summaries of enriched GO terms found in clusters.Blue boxes show GO term summaries related to metabolism, orange boxes related to the immune system, and green boxes related to vitamin D. Other box colors were unique categories.

Figure 5
Figure 5 Differentially methylation (DM) site summary information.(A) Multidimensional scaling (MDS) plots of methylation data(log2(Methylated_counts þ 2) -log2(Unmethylated_counts þ 2)) for all regions or sites in the liver of transgenic and non-transgenic coho salmon reared in either stream or tank environments (expressed in log fold change, logFC).The first MDS plot was generated from all promoter regions, the second was generated from all gene regions, and the last from all methylation sites outside of gene and promoter regions.(B) Venn diagrams of overlapping promoters, genes, or intergenic sites found between different comparisons.Below the promoter and gene Venn diagrams are lines showing the number of overlapping DM regions between promoter and gene based analyses.
R.D. designed the experiment.J.L.L. and K.C. analyzed the data and wrote the paper with L.B., M.C., E.R., and R.D. M.C. and R.D. designed and conducted the sampling.M.C. conducted the qPCR analyses.All co-authors contributed substantially to revisions.

Table 1
Number of differentially expressed genes