Strong Parallel Differential Gene Expression Induced by Hatchery Rearing Weakly Associated with Methylation Signals in Adult Coho Salmon (O. kisutch).

Human activities and resource exploitation led to a massive decline of wild salmonid populations, consequently, numerous conservation programs have been developed to supplement wild populations. However, many studies documented reduced fitness of hatchery-born relative to wild fish. Here, by using both RNA sequencing and Whole Genome Bisulfite Sequencing of hatchery and wild-born adult Coho salmon (Oncorhynchus kisutch) originating from two previously studied river systems, we show that early-life hatchery-rearing environment-induced significant and parallel gene expression differentiation is maintained until Coho come back to their natal river for reproduction. A total of 3,643 genes differentially expressed and 859 coexpressed genes were downregulated in parallel in hatchery-born fish from both rivers relative to their wild congeners. Among those genes, 26 displayed a significant relationship between gene expression and the median gene body methylation and 669 single CpGs displayed a significant correlation between methylation level and the associated gene expression. The link between methylation and gene expression was weak suggesting that DNA methylation is not the only player in mediating hatchery-related expression differences. Yet, significant gene expression differentiation was observed despite 18 months spent in a common environment (i.e., the sea). Finally, the differentiation is observed in parallel in two different river systems, highlighting the fact that early-life environment may account for at least some of the reduced fitness of the hatchery salmon in the wild. These results illustrate the relevance and importance of considering both epigenome and transcriptome to evaluate the costs and benefits of large-scale supplementation programs.

Massive Coho salmon population declines led to the supplementation of hatchery-born fish into the wild. However, despite the apparent absence of genetic differentiation, hatchery fish still display a lower fitness compared with their wild congeners. Here we highlight, in parallel in two distinct river systems, that early-life hatchery environment induces both significant gene expression and epigenomic (i.e., DNA methylation) differentiation that may account for the reduced fitness of hatchery fish in nature. However, the link between patterns of epigenomic and that of gene expression was weak, which suggests that DNA methylation is not the only player in mediating hatchery-related expression differences. Finally, our results also indicate that the dynamic of epigenomic variation should be considered to further evaluate the costs and benefits of large-scale supplementation with captive bred fish.

Introduction

Worldwide exploitation of wild animals and plants by means of fishing, hunting, and other harvesting activities along with rapid environmental changes is rapidly depleting their abundance (Laikre et al. 2010). Numerous conservation programs have thus been developed toward the goal of rehabilitating declining or threatened species through the release of captive or cultivated breeds (Olney et al. 1994; Waples and Drake 2004; Allendorf et al. 2010; Laikre et al. 2010; Allendorf 2017). However, although translocation or supplementation of captive species can lead to demographic increases in population size, they may also come with negative genetic effects on wild populations (i.e., loss of genetic diversity, local adaptation, or local population structure) (Frankham 2008; Laikre et al. 2010; Allendorf 2017) as well as reduction in fitness (Ford 2002; Araki et al. 2007; Araki et al. 2008). Salmonid species are among the most exploited and farmed fish worldwide due to their major socioeconomic value (Aas et al. 2018). Consequently, numerous wild salmonid populations have undergone dramatic declines leading to the establishment of a broad diversity of management and conservation programs (Noakes et al. 2000; Utter 2004; Irvine and Fukuwaka 2011; Aas et al. 2018).

Despite massive supplementation programs, numerous wild salmonid populations are still declining (Augerot 2005; Araki et al. 2008). Many studies have focused on how genetic differentiation induced by selection during hatchery rearing may impact fitness in the wild (Allendorf et al. 2010; Araki and Schmid 2010; Mäkinen et al. 2015; Christie et al. 2016). Although integrated hatchery management programs aim to minimize the potential genetic differentiation, many studies have documented that hatchery-born salmon display maladaptive traits and reduced fitness in the natural environment compared with their wild counterpart despite negligible to weak genetic differentiation between them (Araki et al. 2007; Chittenden et al. 2010; Christie et al. 2014; Neff et al. 2015; Christie et al. 2016; Le Luyer et al. 2017). For example, it has been shown that a single generation of hatchery rearing leads to the differential expression of 723 genes with wild congeners in the steelhead trout (Oncorhynchus mykiss) (Christie et al. 2016). These differences may impact hatchery-born fish fitness as differentially expressed genes were involved in immunity and metabolism (Christie et al. 2016). As a consequence, it is highly relevant to understand the possible role of non-genetic processes underlying such “maladaptive” phenotypic plasticity inducing fitness reduction.

Indeed, it has been well documented that non-genetic sources of heterogeneity may induce phenotypic plasticity (Long et al. 2017; Capp 2021) with important evolutionary consequences, either in an adaptive or maladaptive manner (Aubin-Horth and Renn 2009). In particular, epigenomic variations such as DNA methylation, histone modification, small RNA sequences, or nucleosome positioning (Richards et al. 2017; Laubach et al. 2018) have the potential to modulate gene expression and thus the phenotypic response to environmental stressors (Verhoeven et al. 2016a). Consequently, epigenomic variation may play an important role in evolution by influencing the effect of natural selection and fitness of natural populations (Pfennig et al. 2010; Ecker et al. 2018; Laubach et al. 2018; Strader et al. 2020).

Major questions remain about how environmentally induced epigenomic modifications influence individual gene expression and in turn its potential role for species evolution (Aller et al. 2018; Ecker et al. 2018; Capp 2021; Christensen et al. 2021). For example, both positive and negative correlations between gene expression and DNA methylation have been observed in human cancer cell suggesting both direct and indirect effects of DNA methylation (Long et al. 2017). Additionally, there is conflicting evidence pertaining to genomic regions in regulating expression with DNA methylation in gene promoters (Lowdon et al. 2016), enhancers (Bogdanović et al. 2016), and first introns (Anastasiadi et al. 2018). Recently, Lioznova et al. (2019) demonstrated that methylation at a single CpG site could be more predictive of the level of gene expression than the average promoter/gene body methylation rate; they defined such sites as “CpGs traffic lights” (CpGs TL) (Lioznova et al. 2019).

For Coho salmon (Oncorhynchus kysutch), early-life environmentally induced epimarks (i.e. DNA methylation) during hatchery rearing have been documented both in juveniles (Le Luyer et al. 2017) and in germ line cell of mature males (Leitwein et al. 2021) suggesting a possible role of epigenomic modifications in the observed lower reproductive success of hatchery-born individuals (Neff et al. 2015). However, a mechanistic link between differential methylation, gene expression, and the resulting phenotypic differences observed between wild and hatchery-born individuals (Neff et al. 2015) remains unclear although it has been recently reported in Atlantic Salmon (Salmo salar) that environmental stressors (i.e. increased temperature and hypoxia) induced epigenetic marks that correlated with gene expression, implying epigenetically mediated physiological acclimation to environmental changes (Beemelmanns et al. 2021). This study thus strengthens the hypothesis that hatchery-induced early-life epigenomic marks may affect gene expression and thus be partly responsible for the observed reduced fitness in hatchery-born individuals. However, DNA methylation can be reversible (Edwards et al. 2017) and the dynamics of methylation is still unclear. Consequently, additional research is needed to enlighten how early-life epigenomic modifications may modulate fitness through gene regulation and its potential role in evolution (Angers et al. 2010; Verhoeven et al. 2016; Richards et al. 2017).

Our previous studies on DNA methylation in this system were performed on both white muscle tissue and sperm (Le Luyer et al. 2017; Leitwein et al. 2021). Because it may be tissue specific (Blake et al. 2020; Laporte et al. 2019), documenting patterns of DNA methylation in various types of tissues are necessary toward a full understanding of epigenetic variation in adaptation and evolution (Pfennig et al. 2010; Ecker et al. 2018; Laporte et al. 2019; Strader et al. 2020; Anastasiadi et al. 2021). In this study, we performed tissue-matched RNA sequencing and Whole Genome Bisulfite Sequencing (WGBS) on liver tissue from the hatchery and wild-born adult Coho salmon originating from two previously studied river systems. In doing so, our goals were to assess whether 1) there are significant and parallel differences in gene expression between hatchery and wild-born adult individuals from two different river systems and whether 2) patterns of gene expression were associated with variation in DNA methylation. To do so, we first searched for differentially expressed genes between hatchery and wild-born adult individuals. Then, we performed a gene network analysis to identify networks of coexpressed genes associated with the rearing environment (i.e. wild or hatchery) and their associated biological processes. To investigate the relationship between gene expression and methylation, we assessed methylation at two levels; first, we estimated the median promoter/gene body methylation. Then, we used single dinucleotide CpG methylation within the associated promoter/gene body to search for CpG TL (i.e. CpGs with a significant correlation between methylation and the associated gene expression) (Lioznova et al. 2019). Finally, we compared genes with a significant correlation between methylation and gene expression for both methods and performed a gene ontology analysis to retrieve the biological processes being involved.

Results

Differentially Expressed Genes

RNA sequencing performed on the 47 Coho salmon livers generated an average of 57.4 million pair-end reads per sample (range: 33.6–118.7 M) and expression was quantified for 25,246 genes that passed minimum coverage requirements. We used a negative-binomial generalized linear model as implemented in edgeR (Robinson et al. 2010) to test for differential expression due to factors of rearing environment, sex, river, and all possible interactions. Overall, no gene showed significant interactions between sex and rearing environment or between sex and river, whereas, eight genes had a significant interaction between river and rearing environment. Using global false discovery rates (FDRs), we identified 3,643 genes that were differentially expressed in parallel between hatchery and wild-born Coho from Conuma and Quinsam river systems (figs. 1 and 2). A larger proportion of genes (χ2df = 1 =9.2, P-value = 0.002) were downregulated (52.5%; 1,913 genes) in hatchery fish compared with wild fish (47.5%; 1,730 genes). We also identified 4,801 genes that were differentially expressed between populations from the two rivers (fig. 1). A larger proportion of genes (χ2df = 1 = 8.7, P-value = 0.003) were upregulated (52.1%; 2,503 genes) in the Conuma River population compared with that of the Quinsam River population (47.9%; 2,298 genes). Sex had the strongest main effect on gene expression with 8,852 genes differentially expressed between males and females. Overexpression of genes (51.3%; 4,537 genes) was significantly more common in males (χ2df = 1 = 5.6, P-value = 0.03) than females (48.7%; 4,315 genes). There was considerable overlap of genes (N = 3,943) displaying the influence of multiple main effects on their expression (fig. 1). This degree of overlap was statistically much higher than would be expected by chance (N permutations: 1,000, mean ± SD overlap expected: 3,166 ± 34.9, P < 0.001). This pattern was primarily driven by the excessive overlap of differentially expressed genes (N = 1,444; fig. 1) identified between rivers and rearing environments (N permutations: 1,000, mean ± SD overlap expected: 692 ± 21.8, P < 0.001).

Fig. 1.

Patterns of shared differentially expressed (DE) genes detected across the three main effects of rearing environment, sex, and river.

Fig. 2.

Patterns of DE genes (N = 3,643) in the liver between hatchery- and wild-reared Coho salmon. Expression values are TMM-normalized and log-transformed CPM sequencing reads. Samples and genes are clustered on the x- and y-axes, respectively, using Euclidean distance.

Gene Network Analysis

Weighted gene coexpression network analysis (WGCNA) identified 25 modules of coexpressed genes in liver tissue. Twenty coexpressed modules exhibited correlations (P < 0.05) with at least one of the explanatory variables (i.e. rearing environment, sex, river; supplementary fig. S1, Supplementary Material online). Among these, five modules showed a significant correlation with the rearing environment (FDR < 0.05; fig. 3). The module having the strongest correlation with rearing environment (i.e. darkorange, r = −0.61, P = 5.8 × 10−6; fig. 3) comprised 859 genes and was enriched for the immune system (i.e. tumor necrosis factor [TNF] and Interleukin [IL]-1β signaling and inflammatory responses) as well as circadian rhythm and response to reproductive hormone (i.e. estradiol and progesterone) biological processes (supplementary table S1, Supplementary Material online). The five modules all showed a negative correlation with the rearing environment (fig. 3), which indicated that globally, genes of this module are downregulated in hatchery-born adult Coho compared with adult wild-born Coho. For the darkorange module, there was a strong correlation between gene module membership (correlation with the main axis of module variation) and gene significance (absolute value of the correlation between hatchery origin and gene expression) for genes in the dark orange module (r = 0.63, P < 0.001, fig. 3). The darkorange module contained 443 of the genes (51.6% of the module) that were also identified as being differentially expressed between rearing environments. The other four modules were enriched for biological functions in phagocytosis (coral), endosome to lysosome transition (lightcoral), regulation of growth and differentiation of skeletal muscle (bisque4), and mitochondrial gene expression as well as overlap with elements of the immune signaling pathways identified in the darkorange module (lightsteelblue1) (supplementary table S1, Supplementary Material online).

Fig. 3.

WGCNA; (A) representation of the five modules of the coexpressed genes in liver tissue and their correlation coefficient with the rearing environment and P-values in brackets. (B) Representation of the 859 genes of the darkorange module membership (correlation with the main axis of module variation) and the gene significance which is the absolute value of the correlation between hatchery origin and gene expression.

Methylation

WGBS was performed for 47 individuals and produced a mean of 225 million of reads per individual. A mean of 147 million of reads per individual was mapped to the reference genome masked for cytosine to thymine (C–T) polymorphism. Finally, a mean of 17 million symmetric CpG context cytosines with at least 10× coverages were retrieved per individual.

Relationship Between Methylation and Gene Expression

Median Promoter/Gene Body Methylation

Linear mixed models (LMMs) were built to evaluate the relationship between gene expression of the 3,643 differentially expressed genes and the 859 genes belonging to the darkorange module and median methylation across gene body (hereafter: median gene methylation) and putative promoter regions (5 kb region before the gene body; hereafter: median promoter methylation). For the 3,643 genes identified as differentially expressed between hatchery and wild fish, we found 10 genes for which gene expression was associated with median methylation within gene body (table 1, fig. 4, supplementary fig. 2, Supplementary Material online). Seven and three genes displayed a negative and a positive correlation between the gene expression and the median gene methylation, respectively, with a slope coefficient ranging from −0.354 to 0.631 (table 1 and fig. 4). We also detected three genes for which gene expression was positively correlated with median methylation within the putative promoter region, with a slope coefficient ranging from 0.727 to 1.992 (table 1, fig. 4). For the 859 genes comprised within the darkorange module, 16 genes displayed a significant correlation between gene expression and methylation within the gene body (table 1, fig. 4, supplementary fig. 6, Supplementary Material online). Only two of these genes were identified in the differential expression analysis. Eleven and five genes displayed a positive and a negative correlation, respectively, with a slope coefficient ranging from −1.024 to 0.515 (table 1, fig. 4).

Table 1

Summary of the relationship between gene expression and methylation of the differentially expressed (DE) genes between wild and hatchery-born individuals and the genes from the darkorange module (DO)

Analyses	Region	Gene_id	Intercept	Slope	P-value	R 2m	FDR	Mean meth	Sd meth	Mean logCPM	Sd logCPM	N CpG TL
DE and DO	GeneBody	LOC109865485	4.801	−0.354	4.43E−05	0.149	0.040	0.231	0.210	6.808	0.968	0
DE	GeneBody	LOC109867584	0.769	−0.299	2.07E−04	0.204	0.081	0.766	0.214	1.099	0.900	2
DE	GeneBody	nfkbie	−0.239	−0.299	5.70E−05	0.198	0.041	0.531	0.387	−0.353	0.897	0
DE	GeneBody	LOC109898445	2.122	−0.287	6.60E−06	0.265	0.024	0.830	0.156	3.053	0.720	2
DE	GeneBody	LOC109902736	0.880	−0.283	1.43E−05	0.231	0.026	0.721	0.176	1.308	0.759	1
DE	GeneBody	LOC109864871	1.302	−0.243	2.57E−05	0.129	0.031	0.806	0.227	1.858	0.758	3
DE	GeneBody	LOC109893103	0.774	−0.242	2.26E−04	0.172	0.081	0.734	0.218	1.055	0.709	0
DE and DO	GeneBody	LOC109898379	2.917	0.138	1.01E−04	0.129	0.052	0.771	0.327	4.241	0.454	0
DE	GeneBody	LOC109874102	−0.658	0.250	7.26E−05	0.187	0.043	0.775	0.330	−0.977	0.724	1
DE	GeneBody	LOC109895975	2.424	0.631	1.61E−04	0.157	0.072	0.564	0.428	3.861	2.044	0
DE	Promoter	LOC109870111	1.992	0.136	8.64E−05	0.147	0.100	0.637	0.334	2.856	0.437	0
DE	Promoter	nek9	1.886	0.187	2.89E−05	0.175	0.050	0.564	0.418	2.707	0.500	0
DE	Promoter	LOC109870513	0.727	0.216	1.80E−05	0.247	0.050	0.666	0.359	1.072	0.569	0
DO	GeneBody	LOC109865878	4.614	−0.215	6.67E−04	0.121	0.056	0.076	0.121	6.653	0.718	0
DO	GeneBody	LOC109903482	−1.818	−1.024	5.83E−13	0.685	0.000	0.674	0.267	−2.567	1.680	6
DO	GeneBody	LOC109887160	2.675	−0.471	1.77E−03	0.115	0.099	0.859	0.103	3.792	1.767	0
DO	GeneBody	LOC109885345	2.506	−0.263	1.87E−03	0.145	0.099	0.626	0.360	3.568	0.938	0
DO	GeneBody	LOC109907184	2.927	0.091	1.28E−03	0.164	0.090	0.836	0.278	4.233	0.291	1
DO	GeneBody	LOC109882707	1.039	0.110	5.61E−04	0.218	0.053	0.585	0.413	1.501	0.307	1
DO	GeneBody	jak1	2.668	0.137	1.63E−03	0.124	0.099	0.772	0.312	3.867	0.485	0
DO	GeneBody	dhx57	2.041	0.185	1.57E−03	0.124	0.099	0.936	0.080	2.954	0.686	0
DO	GeneBody	LOC109889095	2.040	0.185	5.08E−04	0.219	0.053	0.553	0.432	2.996	0.564	1
DO	GeneBody	LOC109909164	1.085	0.188	7.27E−04	0.125	0.056	0.816	0.286	1.558	0.598	0
DO	GeneBody	LOC109874412	1.494	0.306	4.99E−04	0.111	0.053	0.807	0.242	2.174	1.156	1
DO	GeneBody	LOC109871914	1.376	0.332	1.85E−04	0.166	0.026	0.810	0.291	2.018	1.000	0
DO	GeneBody	LOC109906505	2.086	0.438	1.17E−04	0.201	0.020	0.820	0.209	2.961	1.255	11
DO	GeneBody	LOC109883217	−0.928	0.515	5.36E−05	0.256	0.015	0.803	0.209	−1.346	1.415	1

Region, region of the analyses (i.e., gene body or promotor); Gene_id, gene identification; Intercept, value of the model intercept; Slope, value of the slope coefficient; R²m, value of the marginal correlation; FDR, false discovery rate; Mean meth, mean methylation of the region; Sd meth, standard error of the mean methylation; Mean logCPM, mean of the gene expression (in count per million); Sd logCPM, standard error of the mean gene expression; N CpG TL, number of CpGs TL detected.

Fig. 4.

Relationship between methylation and gene expression of the 3,643 DE genes between wild and hatchery-born individuals (A) slope coefficient of the generalized linear model (GLM) of the median gene methylation (orange circle) and median gene promotor methylation (5 kb; blue circle) and the gene expression models, (B) slope coefficient of the GLM of CpG Traffic Light methylation and the expression of the associated genes (orange circle) and the associated promotor regions (blue circle). (C) Histogram of the number of CpGs TL per DE genes. Relationship between methylation and gene expression of the 859 coexpressed genes of the darkorange (DO) module: (D) slope coefficient of the GLM models of the median gene methylation (orange circle) and the gene expression, (E) slope coefficient of the GLM models of the CpG Traffic Light methylation and the expression of the associated genes (orange circle) and the associated promotor regions (blue circle). (F) Histogram of the number of CpGs TL per DO genes.

CpG “Traffic Light”

A total of 316,244 filtered CpGs within the gene body of the 3,643 differentially expressed genes and 58,847 filtered CpGs within their putative promoter regions were kept for further analysis. For the 859 genes within the darkorange module, a total of 69,365 and 13,501 filtered CpGs within the gene body and their putative promoter regions, respectively, were kept for further analysis. As for the average methylation analysis, we used LMMs to evaluate the relationship between methylation l and gene expression levels. CpGs for which methylation significantly explained gene expression (FDR < 0.1) were considered as CpGs TL. A total of 527 and 12 CpG TL were found within the gene body and the promoter region of 399 and 11 genes differentially expressed between wild and hatchery with a slope coefficient ranging from −0.80 to 0.61 and −0.76 to 0.51, respectively (fig. 4). A total of 311 genes contained 1 CpG TL, 61 genes contained 2 CpGs TL, 22 genes contained 3 CpG TL, 2 genes contained 4 CpG TL, and 3 genes contained 5, 7, and 9 CpGs TL (fig. 4). In the putative promoter region, we found nine genes with one CpG TL and two genes with three CpG TLs within the 5 kb before the gene body.

A total of 142 CpG TLs were observed among 95 genes of the darkorange module and 1 within the promotor region with a slope coefficient ranging from −1.05 to 0.62 and 1.47, respectively (fig. 4). A total of 74 genes contained 1 CpG TL, 15 genes contained 2 CpGs TL, 1 gene with either 3, 5, 7, and 11 CpG TLs, and 2 genes with 6 CpG TLs (fig. 4).

Genes Found by Both Methods

Five of the 10 genes displaying a significant correlation between gene expression and the median methylation of the gene body also contained at least one CpG TL (table 1). The median CpG methylation, gene expression, gene models, and CpG TL models of these five genes are presented in Supplementary figs. 3–5, Supplementary Material online (see table 1 for gene names). For example, gene LOC109864871 on chromosome 19 displayed the strongest R² (R²c = 0.759, table 1) and contained three CpG TLs (table 1). The distribution of the CpGTL and the 10 genes are represented in fig. 5.

Fig. 5.

Slope coefficient of the GLM models showing the relationship between gene expression of DE genes (violet) and DO modules genes (orange) and methylation of each CpG TL (diamond) and median gene body methylation (square). The relative position of each CpG TL (diamond) and genes (square) is reported for each 30 Coho salmon linkage groups.

For the darkorange genes analysis, we found 7 of the 16 genes that displayed a significant correlation between the gene expression and the median methylation of the gene body also contained at least one CpG TL (table 1). Median methylation, gene expression, gene models, and CpG TL models of those seven genes are presented in supplementary figs. 6–9, Supplementary Material online (please refer to table 1 for gene names).

Gene Ontology Terms

CpG TLs located in gene bodies were associated with genes enriched for biological functions in brain and nervous system development, regulation of synaptic transmission, and response to stimulus (supplementary table S2, Supplementary Material online). CpG TLs located in gene promoter regions were enriched for genes associated with response to testosterone and long-term memory. In contrast, median gene methylation that correlated with gene expression was associated with genes enriched for immune system processes including regulation of interferon gamma production (supplementary table S3, Supplementary Material online). Median promotor methylation that correlated with gene expression was associated with genes enriched for telomerase-related functions and hyperosmotic response (supplementary table S3, Supplementary Material online).

Discussion

Dramatic declines of wild salmonid populations in the 1990s led to the establishment of conservation and hatcheries programs with the goal of restoring depleted populations. Despite these efforts, and hatchery programs specifically managed to minimize genetic differences between wild and captive bred fish, hatchery-born individuals generally display a lower reproductive success than their wild counterparts in natural habitats (Christie et al. 2014; Neff et al. 2015). Lower relative fitness of hatchery fish in the natural environment (Araki et al. 2008) has been reported even after a single generation of hatchery rearing (Araki et al. 2007; Christie et al. 2016). Here, we observed pronounced differential gene expression in parallel between hatchery and wild-born adult Coho salmon from two river systems after almost 2 years spent growing in the same environment (i.e., at sea) and about 1.5 years spent in the hatchery environment. The extent of parallel differential gene expression between wild and hatchery Coho was on the same order (80%, i.e., 3,643 vs. 4,537 differentially expressed genes) than that observed between Coho from these two river systems. Moreover, parallel downregulation of gene expression in hatchery-born individuals occurred despite no detectable genome-wide genetic differentiation between hatchery and wild fish from a given river (Le Luyer et al. 2017; Leitwein et al. 2021). Differential expression targeted important biological processes related to the immune response that were downregulated in hatchery-born fish, thus strengthening the hypothesis that early-life environment may have a significant impact on individual fitness. It has been documented that variation in the level of gene expression variation can be attributed to epigenomic variations (Elowitz et al. 2002; Pujadas and Feinberg 2012; Ecker et al. 2018) underlying genotype-by-environment interactions (Christensen et al. 2021). Yet, we did not detect broad associations between differential DNA methylation and the observed gene expression differences in this study based on WGBS. Nevertheless, we observed several significant relationships between gene expression and the whole gene body methylation as well as methylation at a single CpGs. These results strengthen the emerging view that the dynamic control of gene expression by DNA methylation is not straightforward and involves complex, rather than direct associations (Jones 2012; Baubec and Schübeler 2014; Ambrosi et al. 2017). Below we discuss the relevance of our results for understanding the fitness consequences of hatchery rearing and the role that methylation plays in mediating the underlying gene expression differences.

Early-Life Environment-Induced Differential Gene Expression

The hypothesis of a negative hatchery-rearing effect on individual fitness has previously been supported by the observation of pronounced gene expression differences between recently hatched hatchery and wild juvenile Rainbow Trout (O. mykiss) in the absence of genetic differences between them (Christie et al. 2016). Here, we found 3,643 genes that were differentially expressed in parallel between wild and hatchery-born adults from two different river systems. This is particularly noteworthy as the fish in our study first experienced hatchery conditions during the first year and half of their existence, followed by more than 18 months in a common growing environment (i.e., the sea). This result provides compelling evidence that environmental exposure during early developmental stages exerts a strong and lasting effect on the physiology of individuals and likely has a non-negligible influence on their fitness (Le Luyer et al. 2017; Leitwein et al. 2021).

The gene coexpression module that had the strongest association with rearing environment (i.e., darkorange: 859 genes) was enriched for immune systems response processes (specifically early inflammatory responses mediated by TNF and IL-1 signaling) that were downregulated in hatchery-born individuals. Reduced disease resistance for hatchery-reared Coho salmon has previously been reported (Salonius and Iwama 1993) and is consistent with our gene expression results. Differential regulation of genes involved in the immune system has also been reported in hatchery-reared Rainbow Trout (Christie et al. 2016). In Sockeye Salmon (Oncorhynchus nerka), differences in expression of genes involved in the immune system have been linked with the ability to successfully migrate and reproduce, thus potentially impacting fitness (Miller et al. 2011). The consistency of differential expression of immune system processes between hatchery and wild fish across species and life stages, as well as the obvious link between immune system functions and salmon fitness, suggests that hatchery-induced reductions in fitness may be at least partly mediated through a reduction of immune system competency.

Hatchery fish also showed a lower expression of genes involved in the circadian rhythm which could have important fitness consequences. Circadian rhythm may be subject to strong selective environmental pressures as, for example, feeding during day light highly increases the risk of predation and inversely feeding at night may decrease the food availability (Yerushalmi and Green 2009). Alteration of the circadian rhythmicity occurring when food is continuously available in the hatchery environment can lead to maladaptative activity behavior of hatchery fish into the wild (Alioravainen et al. 2019). Finally, differential expression of genes implicated in hormonal production may also directly be related to the individual fitness as a seasonal production of gametes that mismatch the natural hatching period will result in a lower reproductive success (Fregeneda-Grandes et al. 2013; Christie et al. 2014). Altogether, such differential gene expression observed between hatchery and wild fish in different salmonid species and in the absence of apparent genetic differences bring increasing support to the hypothesis that early-life rearing environment could have a significant impact on individual fitness.

Sex and Population Effects on Patterns of Gene Expression

In addition to the differential expression observed due to rearing environments, we reported strong gene expression differences between sexes (8,852 genes) and rivers (4,537 genes). Strong patterns of differential gene expression between sexes were expected as sex-biased gene expression is a widely accepted mechanism for resolving sexual conflict in highly sexually dimorphic species (Mank 2017), which is consistent with the sexual dimorphism and pronounced differences in gene expression observed in salmonids, for instance in Brook Charr, Salvelinus fontinalis (Sutherland et al. 2019). More surprising was the fact that a significant number of differentially expressed genes observed between populations from the two rivers (1,444 genes) was also differentially expressed between rearing environments. The lack of widespread evidence for an interaction effect between river and rearing environment suggests that these genes are differentially expressed in the same direction but to a higher degree in one population over the other. Several processes involved in the response to a viral pathogen were among the most overrepresented among the differentially regulated genes between rivers. A disease outbreak in one population would explain both the river-specific differential expression patterns and the significant overlap between main effects. As discussed above, immune response expression pathways appear to be impaired in hatchery-reared fish and a disease outbreak in one river would generate both river-specific patterns of expression, whereas, altering the nature of expression differences between hatchery and wild fish.

Weak Link Between Cis-Acting DNA Methylation and Gene Regulation

The role of methylation in dynamically regulating gene expression remains debated in the literature, especially because the molecular mechanisms driving expression are not fully understood (Long et al. 2017; Anastasiadi et al. 2018). Here, our results reveal at least some degree of linkage between gene expression and both the gene body median methylation and methylation at single CpGs. Here and in previous studies, by comparing two river systems, we observed significant differences in patterns of methylation between genetically similar hatchery and wild-born Coho (Le Luyer et al. 2017; Leitwein et al. 2021), which indicates that early-life environment has a significant impact on the methylation pattern and may play an important role in evolution by impacting population fitness (Pfennig et al. 2010; Ecker et al. 2018; Laubach et al. 2018; Strader et al. 2020). However, we detected only a small proportion of genes where environmentally induced DNA methylation difference was significantly associated with that of gene expression. These results are broadly similar to those of other recent work in this species that also found weak linkages between methylation and gene expression underlying genotype-by-environment interactions in this species (Christensen et al. 2021). The small proportion of differentially expressed genes for which we observed a correlation between DNA methylation and gene expression suggests that cis-acting DNA methylation is not the primary mechanism mediating the entire breadth of hatchery-related gene expression differences. Whereas, WGBS is one of the more robust methods to detect genome-wide DNA methylation, it is biased toward characterizing high-density CpG regions (Beck et al. 2021) and thus power to detect DNA methylation associations with gene expression in genomic regions of lower CpG density may have been lower than expected.

Moreover, although we reported genes for which hyper-methylated regions are associated with downregulation genes, as classically observed (Edwards et al. 2017), we also observed hyper-methylated regions associated with upregulated genes. Such a relationship has previously been observed in Atlantic Salmon where hypo-methylated promoter CpG sites were positively correlated with downregulation of the associated gene (ucp2) (Beemelmanns et al. 2021). Clearly then, the relationship between methylation and gene expression is not straightforward and more complex than previously believed (Jones 2012; Baubec and Schübeler 2014; Ambrosi et al. 2017). As in Lioznova et al. (2019), we also reported that methylation at a single CpG can be significantly correlated with the associated gene expression and that significant relationship may be missed by only considering whole gene body methylation. Whereas, this indicates that methylation at single CpG sites should be considered in future studies, this also strengthens the fact that patterns of DNA methylation are dynamic and context-dependent (Beemelmanns et al. 2021). We did not find support for methylation-influenced differentially expressed genes to be “core” modules genes (i.e., high intramodule connectivity) with expected pleiotropic effects on the expression of many other genes. This suggests that methylation effects at key genes with subsequent cascading effects on gene expression are not widespread and does not explain the lack of widespread methylation-expression associations.

It is possible that methylation of genomic features outside the regions we characterized have an influence on the observed expression patterns we report. The limited availability of other epigenomic-related resources for Coho salmon (e.g., chromatin accessibility, histone modifications, transcription-factor binding positions, etc.) prevented us from identifying other genomic features (e.g., enhancers) that are known to interact with methylation patterns to assess more complex associations of DNA methylation with gene expression. Indeed, in humans, methylation at enhancers has been reported to correlate better with expression than methylation of promotors (Fleischer et al. 2017). Thus, our associations of DNA methylation and gene expression should be considered a conservative estimate of the influence of DNA methylation on gene expression. Another important consideration is the possibility that hatchery-associated differential DNA methylation influences developmental checkpoints and has long-lasting effects on expression patterns that are not captured in our study. Future work will need to refine our understanding of the sequence and timing of hatchery-related differential methylation as well as cell- and tissue-specific patterns of differential methylation to uncover the specific ways in which the hatchery environment alters fish phenotype.

Implication for Conservation

To conclude, we demonstrate in two Coho populations from distinct river systems, that early-life hatchery environment-induced pronounced and parallel gene expression differentiation (80% of the level observed between rivers) as well as epimarks that persist after 10 months in the same marine growing environment. These results further support the hypothesis that despite hatchery effort to diminish hatchery and wild-born genetic differentiation, for instance, by implementing integrated hatchery program (Araki et al. 2008), early-life hatchery-rearing environment reduced fitness in the wild. Our study also highlights the fact that hatcheries should consider the development of environmental rearing conditions that mimic, as close as possible, wild conditions to limit gene expression differentiation occurring during early developmental stages. On a broader scale, whereas, the link between methylation and gene expression remains weak, dynamics of epigenomic patterns should be considered in the management and conservation of salmonids strategies to evaluate the cost and benefits of large-scale supplementation with captive bred fish.

Methods

Sampling

Genetically, indistinguishable hatchery and wild Coho salmon (Le Luyer et al. 2017; Leitwein et al. 2021) were sampled in the Quinsam and Conuma river systems in British Columbia (Canada). Both Quinsam and Conuma hatcheries are part of Fisheries and Oceans Canada’s Salmon Enhancement Program that has a primary production strategy (see Le Luyer et al. 2017 for more details). Briefly, the aim of these “integrated” hatchery programs is to use all local returning fish, both wild and hatchery fish, as broodstock to minimize genetic differentiation between the two systems. Returning adults were sampled before spawning at Quinsam and Conuma hatcheries, on November 1 and 2, 2017, respectively. For each sampling site, hatchery-born fish were identified by their clipped adipose fine. The livers of 48 adult Coho salmon returning to their natal river after spending about 1.5 years in the ocean were sampled, including 24 fish from Quinsam hatchery (seven wild and five hatchery males, six wild and six hatchery females) and 24 fish from Conuma hatchery (six wild and six hatchery males, seven wild and five hatchery females). Fish were kept in a large holding tank with flowing water before being euthanized and sampled. All samples were stored in RNAlater in the fridge (4 °C) for 24 h and then transferred into a −80°C freezer. Livers were used because of their relative homogeneity compared with other organs and it also displays the strongest gene expression activity (in terms of the number of genes being expressed in this organ) compared with brain, gills, and kidneys in salmonid migratory species (Giger et al. 2007).

DNA Extraction and WGBS

Genomic DNA was extracted following a protocol of universal and rapid salt extraction (Aljanabi and Martinez 1997). DNA quality control, WGBS library preparation, and 100 bp paired-end sequencing on an Illumina HiSeqX (two individuals per lane) were performed at the McGill University and Génome Québec Innovation Centre (Montréal, Ontario, Canada). One individual (one hatchery-born female from Quinsam) did not pass DNA quality control and was excluded from this study, leaving a total of 47 Coho salmon samples for statistical analyses.

RNA Extraction and RNA Sequencing

RNA was extracted from the liver using RNeasy Kits (Qiagen Inc., Toronto, Ontario, Canada) following the manufacturer’s instructions using on the column DNase treatment. RNA quality was verified using Agilent RNA 6000 Nano chips on the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, California, USA). Only samples with an RNA integrity number ≥ 7 were used for further analysis. RNAseq libraries were prepared at the Centre d’Expertise et de Services Génome Québec (Montréal, Québec, Canada) using the NEBNext Ultra II Directional mRNA Library Prep Kit (New England Biolabs Ltd., Whitby, Ontario, Canada). Individually barcoded libraries were pooled and sequenced with 100 bp paired-end sequencing in a single lane of an S4 flow cell on a NovaSeq 6000 (Illumina Inc., San Diego, California, USA).

Differential Gene Expression

RNA sequencing reads were adaptor- and quality-trimmed using fastp v0.20.0 (Chen et al. 2018) requiring a minimum length of 80 bp and implementing the “–poly_X” option to additionally trim trailing polyA homopolymers. Transcript expression levels were quantified with Salmon v1.1.0 using decoy-aware selective mapping (Patro et al. 2017) based on the gene annotation of the Coho salmon reference genome (Okis_V1, NCBI RefSeq: GCF_002021735.1). Transcript quantification data were processed using the tximport v1.12.3 package (Soneson et al. 2015) within the R v3.6.3 statistical environment (Team 2015) to generate estimates of gene expression level, whereas, correcting for sample-specific transcript length differences. We tested for differential expression using negative-binomial generalized linear models with quasi-likelihood tests as implemented in the edgeR v3.26.8 package (Robinson et al. 2010). Models were fit with the same design as the methylation data with main factors of river, sex, and rearing environment, as well as all possible interactions. FDRs (Benjamini and Hochberg 1995) were calculated and differential expression was determined at FDR < 0.05.

We defined coexpressed networks of genes using WGCNA as implemented in the WGCNA v1.69 package in R (Langfelder and Horvath 2008). We used trimmed mean of M values (TMM)-normalized, log-transformed read counts per million (CPM) for all 25,246 expressed unigenes in the calculation of gene–gene adjacencies (correlations). A signed-hybrid adjacency matrix for all pairs of genes was constructed using the robust biweight midcorrelation raised to the power of nine to approximate a scale-free network. The adjacency matrix was then transformed into a topological overlap dissimilarity matrix and a combination of hierarchical clustering and a dynamic tree-cutting algorithm were used to first define and then merge coexpressed modules of genes (minimum module size of 30 genes).

Gene expression modules (i.e. first module eigenvector: first axis of a principal component analysis conducted on the expression of all module genes) were correlated with factors (sex, river of origin, and rearing conditions) using the biweight midcorrelation. Modules with a significant correlation (P < 0.05) were retained for further analysis. For each significant module-trait correlation, we assessed the biweight midcorrelation of each module gene with the phenotype to infer the “gene significance.” Gene significance was then correlated with module membership (biweight midcorrelation of gene expression versus the module eigenvector) to validate module correlations with each factor and identify core genes in the module.

Methylation Calling

The WGBS reads were trimmed and quality filtered (min quality = 25, min length = 100 bp) with fastp (https://github.com/OpenGene/fastp [Chen et al. 2018]). To avoid interpreting false epigenetic variation with existing C–T polymorphism, we masked the reference genome (NCBI assembly GCA_002021735.1; Okis_V1) from C–T polymorphism identified with whole genome resequencing of 20 Coho salmon (940,406SNPs, maf = 0.05) from four British Columbia rivers, using BEDtools maskfasta v2.26.0 (Quinlan and Hall 2010) as in Le Luyer et al. (2017). WGBS trimmed reads were mapped against the masked coho genome with WALT v1.0 (Chen et al. 2016; https://github.com/smithlabcode/walt) by using default parameters and a maximum allowed mapping for a read (-k) of 10. Individual’s symmetric CpG methylation levels were estimated with MethPipe v.3.4.3 (https://github.com/smithlabcode/methpipe). All symmetric CpG sites with less than 10X coverages were removed.

Relationship Between Methylation and Expression

To assess the relationship between gene expression and methylation levels, we first computed for each individual the median of methylation level within the gene body and within the 5 kb region (i.e. putative promoter regions) upstream (5′-end) of the gene body for the 3,643 differentially expressed genes between hatchery and wild-born Coho and the 859 genes within the darkorange module which was the module having the strongest correlation with rearing environment (see “Results” section). We then built LMMs where the TMM-normalized log-transformed CPM gene expression was introduced as the response variable and the standardized (i.e., center-reduced) median methylation was included as the explanatory variable. River of origin, rearing environment, and sex were incorporated as random effects in the model. P-values were computed to evaluate the significance of the models. Correction for multiple testing was corrected by performing the Benjamini and Hochberg FDR correction (FDR < 0.1) (Benjamini and Hochberg 1995). All analyses were performed using the R package “LME4” (Team 2015).

Fine Scale CpG TL

The link between methylation and gene expression is not straightforward and remains debated (Lioznova et al. 2019). Considering only the median gene and promotor methylation might result in missed local (single CpG) effects of methylation where the signal can be swamped when averaging across the gene body. Following the approach used in a previous study (Lioznova et al. 2019), we examined the relationship between single CpG methylation levels and gene expression by performing LMM for each CpG within the gene body and within the 5 kb putative promoter region before the gene body of the 3,643 differentially expressed genes between wild and hatchery-born Coho and for the 859 genes from the darkorange module. Only CpGs present in at least 90% of the individuals and for which at least 80% of the individuals displayed a methylation value greater than zero were kept to build the models. As previously, the TMM-normalized gene expression in CPM was log + squared transformed and introduced as the response variable, the methylation (standardized) was introduced as the explanatory variable and the river of origin, rearing environment, and sex as random effects. FDR (Benjamini and Hochberg 1995) was applied to P-values of each model. Only CpGs with FDR lower than 0.1 were considered as CpG TL.

GO Terms Analysis

To test for enrichment of biological functions, in particular groups of genes (e.g., differentially expressed, influenced by methylation, module membership, etc.), we first generated Gene Ontology (GO) functional annotation for all genes found in the NCBI gene annotation for the Okis_V1 genome assembly used for mapping (NCBI accession: GCA_002021735.1). We aligned protein sequences for all genes against the SWISS-PROT database (accessed December 2019) using the blastp program with default settings from the BLAST+ v2.6.0 command line application (Camacho et al. 2009). We then used the software Blast2GO v5.2.5 (Conesa et al. 2005) to map GO terms from the top 20 blast hits for each sequence requiring a e-value < 1e−5 to retain matches. Using this annotation, we performed enrichment tests with the topGO v2.38.1 package in R (Alexa and Rahnenführer 2021). Enrichment of GO biological functions was tested based on the “weight01” (Alexa et al. 2006) algorithm and Fisher’s Exact Tests. Unless specified, the entire set of expressed genes (N = 25,246) was used as the reference gene set.

Click here for additional data file.