Environment-driven reprogramming of gamete DNA methylation occurs during maturation and is transmitted intergenerationally in Atlantic Salmon

Abstract An epigenetic basis for transgenerational plasticity in animals is widely theorized, but convincing empirical support is limited by taxa-specific differences in the presence and role of epigenetic mechanisms. In teleost fishes, DNA methylation generally does not undergo extensive reprogramming and has been linked with environmentally induced intergenerational effects, but solely in the context of early life environmental differences. Using whole-genome bisulfite sequencing, we demonstrate that differential methylation of sperm occurs in response to captivity during the maturation of Atlantic Salmon (Salmo salar), a species of major economic and conservation significance. We show that adult captive exposure further induces differential methylation in an F1 generation that is associated with fitness-related phenotypic differences. Some genes targeted with differential methylation were consistent with genes differential methylated in other salmonid fishes experiencing early-life hatchery rearing, as well as genes under selection in domesticated species. Our results support a mechanism of transgenerational plasticity mediated by intergenerational inheritance of DNA methylation acquired late in life for salmon. To our knowledge, this is the first-time environmental variation experienced later in life has been directly demonstrated to influence gamete DNA methylation in fish.


Introduction
The inheritance of environmentally induced epigenetic variation (e.g., DNA methylation, chromatin structure, small RNAs) has been proposed as a mechanism facilitating transgenerational plasticity (Richards 2006;Bell and Hellmann 2019). As a chemical modification of nucleotide bases, DNA methylation has a clear mechanism for multigenerational transfer, however, our current understanding of its role in intergenerational (i.e., parentoffspring transfer) or transgenerational (i.e., multi-generational transfer) epigenetic inheritance in animals is hindered by a lack of data from a wider diversity of organisms and thus evidence remains limited and controversial (Heard and Martienssen 2014;Horsthemke 2018;Skvortsova et al. 2018;Anastasiadi et al. 2021). In mammals, DNA methylation undergoes extensive reprogramming (de-methylation and re-methylation) following fertilization and then again when germ-line tissue differentiates from somatic tissue (Heard and Martienssen 2014). This generally limits the potential for intergenerational inheritance as methylated sites must escape two rounds of re-patterning in each generation. Other well-studied traditional animal models (i.e., worms and flies) lack widespread methylation and thus currently provide no relevant insight (Dunwell and Pfeifer 2014;Hu et al. 2015). However, a fish model, specifically zebrafish (Danio rerio), does not exhibit global erasures and reprogramming of methylation neither following fertilization where the maternal pattern is reprogrammed to match the paternal pattern (Jiang et al. 2013;Potok et al. 2013) nor during germline differentiation Skvortsova et al. 2019) suggesting a greater potential for DNA methylation-mediated epigenetic inheritance in this species. While evidence from medaka (Oryzias latipes) suggests it may also undergo methylation reprogramming similar to mammals (Wang and Bhandari 2019), the lack of consistent patterns continues to make teleost fishes an interesting group in which to investigate the potential for intergenerational epigenetic inheritance.
Environmentally induced DNA methylation variation in a wider selection of teleost fishes has been reported as a result of differences in early developmental environments (Artemov et al. 2017;Le Luyer et al. 2017;Metzger and Schulte 2017;Gavery et al. 2018;Nilsson et al. 2021). These signals stably persist until adulthood (Metzger and Schulte 2017;Leitwein et al. 2021), occur in germ cells (Gavery et al. 2018; Rodriguez Barreto et al. 2019), and a growing body of literature demonstrates that intergenerational transmission can occur (Ryu et al. 2018;Rodriguez Barreto et al. 2019;Berbel-Filho et al. 2020;Heckwolf et al. 2020), thus supporting an overall mechanism of DNA methylation-mediated intergenerational epigenetic inheritance for this clade.
The timing for germ-line incorporation of environmentally induced epigenetic variation is historically believed to be limited to early developmental stages as a result of the separation of the germline and soma, or the so-called "Weismann Barrier" (Monaghan and Metcalfe 2019). Once germline epigenetic reprogramming is completed, this barrier limits somatic influence on the germline and thus, would theoretically be expected to prevent environmentally derived epigenetic changes from being incorporated into gametes and passed to future generations. Yet, emerging evidence suggests that the barrier is more permeable than previously thought (Eaton et al. 2015), with evidence for small RNA mediated epigenetic inheritance (Sciamanna et al. 2019;Duempelmann et al. 2020) and at least one example of differential methylation in mice (Dias and Ressler 2014). However, wide taxonomic evidence for these phenomena is lacking and this makes it difficult to assess the potential for intergenerational transmission of environmentally induced DNA methylation variation acquired after early developmental stages.
Salmonid fish hatcheries provide relevant systems in which to test hypotheses regarding intergenerational epigenetic inheritance. Hatcheries have been used for decades to enhance, supplement, or recover salmonid fish populations (Naish et al. 2007), but have negative consequences for the fitness of hatchery-reared fish in the wild (Araki et al. 2008;Christie et al. 2014;O'Sullivan et al. 2020) presumably due to domestication effects. Numerous studies have failed to demonstrate significant genetic difference between hatchery-origin and natural-origin (wild) salmon (Christie et al. 2016;Le Luyer et al. 2017;Gavery et al. 2018) despite documented pronounced differences in gene expression (Christie et al. 2016). In contrast, DNA methylation divergence has been reported between hatchery-origin and wild salmon for several species (Le Luyer et al. 2017;Gavery et al. 2018Gavery et al. , 2019Rodriguez Barreto et al. 2019;Leitwein et al. 2021). Evidence for intergenerational effects has recently emerged (Rodriguez Barreto et al. 2019), although solely in the context of early developmental hatchery exposure where salmon are artificially reproduced and the offspring reared in the hatchery environment for a period of time before release into the natural environment. Assuming the Weismann Barrier exists in salmon, environmentally induced epigenetic differences would be expected to arise shortly following fertilization and before the differentiation of primordial germ cells.
Alternative hatchery rearing techniques for salmon, including juvenile (here: smolt) to adult supplementation (hereafter SAS; as per Fraser 2016) and live gene-banks, are increasingly being applied to conserve and recover the most critically endangered salmon populations (O'Reilly and Doyle 2007;Stark et al. 2014). In these conservation strategies, wild-born juvenile fish are collected from the wild, reared to adulthood in a hatchery environment, and at the onset of maturation they are again released to the wild to spawn naturally. These strategies are believed to reduce the risks to wild populations that are associated with the use of hatchery technology by allowing natal stream imprinting by juveniles as well as homing, natural mate-choice, and spawning site selection by adults. While these approaches show promise for demographic recovery of some populations (Berejikian and Van Doornik 2018), fitness-related differences between SAS and wild salmon have also been documented (Berejikian et al. 2001(Berejikian et al. , 2005. In these contexts, it is unclear whether environmentally derived epigenetic modifications are a plausible explanation for the observed phenotypic effects given that germ-line DNA methylation patterns, assuming an impermeable Weismann Barrier, should have been determined well before fish experienced a captive environment. As such, we currently lack knowledge of the potential for epigenetically mediated intergenerational effects that may result in heritable declines in fitness in these contexts. These systems also provide opportunities to test the potential for intergenerational transmission of late-life, environmentally induced, DNA methylation variation. In this context, the goal of our study was to test the hypothesis that environmental differences experienced by adults during maturation alter DNA methylation of gametes. Because the zebrafish methylome is reprogrammed to the paternal state (Skvortsova et al. 2019), we first compared DNA methylation in sperm cells from maturing male SAS fish to wild male salmon from overlapping cohorts that had spent 1 year (i.e., grilse) to 2 years (i.e., multi-sea-winter salmon) in the marine environment. We further tested the hypothesis that adult rearing environments influence offspring methylation patterns and that these differences influence offspring phenotypes by creating pure-type crosses of SAS and wild adults and rearing the offspring in a common environment. We characterized growth-related phenotypic differences in 10-month-old juveniles as well as DNA methylation patterns in their livers, to determine the presence of inherited DNA methylation patterns and their potential influence on proxies of fitness.

Methods
Wild Atlantic Salmon (Salmo salar) juveniles ( Figure 1A) were produced naturally by their parents and reared in the tributaries of 60 A B C E D Figure 1 Smolt-to-adult supplementation in the Northwest Miramichi River. Natural origin (wild) juvenile salmon (A) are captured during their migration to the ocean (B) and reared in captivity until adulthood (SAS) at the Miramichi Salmon Conservation Centre (C). Wild salmon continue their marine migration spending 1-3 years feeding in the Labrador Sea (D) before returning to the Miramichi River to spawn. Wild salmon from the same cohort as SAS salmon were captured in various headwater pool habitats during their return migration (E).
the Northwest Miramichi River, New Brunswick, Canada for 2 or 3 years. The smolt run is typically composed of approximately equal proportions of 2-and 3-year-old smolts with a small proportion (<5%) of 4-year-old fish (Chaput et al. 2016). We collected juveniles (smolts) that were migrating to the ocean using a rotary screw trap from May to June 2015 from the mainstem of the near the mouth of Trout Brook ( Figure 1B). They were transported to the Miramichi Salmon Conservation Centre (MSCC, South Esk, NB; Figure 1C) where they were held in 300,000 L tanks on fresh well-water under natural photoperiod until 2017 (i.e., smolt-toadult supplementation or SAS). Water temperature ranged from 5 C to 7 C except during the months of May to September where heat-exchangers were used to extract heat from surface water (Stewart Brook) to warm well-water to between 10 C and 15 C tracking the typical summer increase in temperature. SAS fish were initially fed chopped frozen krill for one month and gradually weaned onto standard hatchery pellet food (Skretting Canada, St. Andrews, NB) by coating it with puré ed krill for an additional 3-4 weeks before feeding solely pellet food. Densities increased throughout the study as fish grew but never exceeded 10 kg/m 3 .
Natural origin (wild) adult Atlantic Salmon returning to spawn in the Miramichi River were collected by seining in September of 2017 and 2018 as a part of regular brood-stock collection program by the Miramichi Salmon Association staff from selected pool habitats in the upper reaches of a major branch (the Little Southwest Miramichi River) of the Northwest Miramichi River ( Figure 1E). Adult fish were transferred to MSCC and held in 5000 L tanks for up to three weeks under natural photoperiod in the same water conditions as SAS fish (well-water). Fish were fasted during holding period and densities in these tanks were maintained between 5 and 7 kg/m 3 .
In 2017, both SAS and wild adult fish were haphazardly netted from holding tanks, gametes were collected, and eggs were artificially fertilized to create pure-type breeding crosses (SAS Â SAS and wild Â wild; N ¼ 8 crosses each). All SAS males and females had spent two winters in the hatchery. All wild males and one female were one sea winter fish (1SW) and the remaining females were 2SW fish. Fertilized eggs were incubated in flow-through troughs until first feeding after which offspring from multiple families from each treatment were mixed and transferred to one of three 5000 L square tanks. Fish were fed to satiation on hatchery pellet food until they reached a size of approximately 1.5 g (mid-August 2018, approximately 10 months, maximum density of 0.7 g/L). Egg sizes did not differ between wild and SAS adults (Mean egg diameter: Wild ¼ 6.8 mm, 95% CI: 6.5-7.1; SAS ¼ 6.6 mm, 95% CI: 6.3-6.9; F 1,6 ¼ 0.82, P-value ¼ 0.4) and survival to the eyed-egg stage was over 90% for all but one SAS family that averaged only 47% survival across three replicate egg baskets (D. Roth unpublished data).
Juvenile offspring fish were netted haphazardly from the tanks, euthanized in an overdose solution of eugenol (Sigma-Aldrich Canada, Oakville, ON, Canada), weight and length measurements taken, and liver tissue was dissected and preserved in RNAlater. We chose liver tissue because it is an organ known for its important role in regulating growth and metabolism (Trefts et al. 2017). Liver also represents a tissue in which a very large number of genes are being expressed in salmonids that could be regulated by epigenetic mechanisms (e.g., Rougeux et al. 2019;Sutherland et al. 2019;) as well as its relative homogeneity of celltypes that could confound methylation analyses (Jaffe and Irizarry 2014). Juvenile samples were genotyped with a panel of 188 SNPs (KASP SNP assays; LGC Biosciences, Beverley, MA, USA) and assigned to their family of origin using COLONY v2.0.6.5 (Jones and Wang 2010). In addition, we determined genetic sex of the offspring using a previously described assay (King and Stevens 2019). We chose two genetically male individuals from each of four pure-type families for each treatment group (N ¼ 16) to conduct whole-genome bisulfite sequencing.
Gamete samples were unavailable from the adults used in the breeding crosses in 2017. The following year (i.e., 2018), we collected sperm samples from eight wild and eight SAS males. All wild males were grilse having spent one winter at sea, while all SAS males had spent two winters in the hatchery. Sperm was collected in 2-mL tubes and stored on ice for between 2 and 6 h. Sperm (250 lL) was centrifuged at 7000 rpm for 10 min, the supernatant discarded, and isolated sperm cells preserved in 1.5 mL of RNAlater. While this experimental design does not allow direct inference of parent-offspring transmission of DNA methylation, it does still allow us to infer intergenerational transfer of DNA methylation due to the fact that we controlled the rearing environment for the juveniles (common garden) and thus these fish only differ by their respective parents' environmental exposures. Because the 2018 fish originated from the following cohort as the 2017 parents, this also allows a degree of repeatability of the effects to be tested.
To characterize C-T polymorphisms that could bias methylation estimates, we combined equal proportions of DNA from all individuals (N ¼ 36) and sequenced them as a pool in one lane of an Illumina HiSeqX. Raw sequencing reads were trimmed using fastp v0.19.9 and then aligned to the Atlantic Salmon genome with bwa mem (Li 2013). Duplicate reads were removed using MarkDuplicates and overlapping 3' ends of paired reads were clipped using the clipOverlap function of BamUtil v1.0.14 (Jun et al. 2015). We called SNPs using a frequency-based approach in freebayes v1.3.1 (Garrison and Marth 2012) that required variant sites to be covered by a minimum of 10 reads and have a minimum of two reads supporting the alternate allele. We retained both C-T and A-G (i.e., C-T on the minus strand) polymorphisms and removed these sites from the methylation results using bedtools v2.26.0 (Quinlan and Hall 2010).
We quantified differential methylation at CpG sites covered by at least one read in all samples where we additionally required a minimum of five reads and a maximum of 20 reads (approximately 99.9th quantile) for at least 12 of 16 juveniles or eight of 12 adults. The sequencing performance for four adult samples (i.e., two HiSeqX lanes) representing two SAS and two wild fish was poor and thus these samples were excluded to reduce biases in methylation estimates due to low coverage for these samples. The minimum coverage filter ensured that differences in methylation were not due to spurious differences in coverage between groups and the maximum coverage filter removed highly repetitive regions where the confidence in mapping accuracy is low. Differential methylation of CpG cytosines (DMC) was determined using beta-binomial models implemented in the DSS v2.32.0 package (Feng et al. 2014) in the R statistical environment v3.6.1 (R Core Team 2019). Methylation levels were first smoothed using a window size of 500 bp and models were fit with group-specific dispersion estimates as implemented in DSS. False discovery rates (FDRs) were calculated according to Benjamini and Hochberg (1995). Differentially methylated regions (DMRs) were determined based on attributes of DMCs, where regions were required to be a minimum of 50 bp long, have > 3 CpGs, and greater than 50% of the CpG sites with a P-value < 0.001. Due to the large number of small contigs in the Atlantic Salmon genome (i.e., >230,000; ICSASG v2; NCBI RefSeq: GCF_000233375.1; Lien et al. 2016), we restricted our analyses to the 29 full-length chromosomes and contigs larger than 10 kb in length (>96% of the ungapped length of the genome).
To determine potential functional consequences of methylation differences we used bedtools v2.26.0 (Quinlan and Hall 2010) to identify gene features associated with DMRs. NCBI RefSeq gene annotation information for the salmon genome was retrieved and genes were associated with differential methylation if any DMRs were located within 5000 bp of a coding region consistent with previous work in salmonids (Le Luyer et al. 2017). Gene ontology information for annotation genes was obtained from the Ssal.RefSeq.db v1.3 (https://gitlab.com/cigene/R/Ssa.RefSeq.db; accessed: June 22, 2020) R package. We tested for enrichment of biological functions for genes associated with DMRs using Fisher's Exact Tests and the 'weight01' (Alexa et al. 2006) as implemented in the TopGO v.2.38.1 package (Alexa and Rahnenfuhrer 2019).
We used a network-based approach to investigate associations of correlated methylation signatures with juvenile phenotypes (i.e., length, weight, and condition factor). We first summarized methylation in nonoverlapping 100 bp windows across the genome. Windows required a minimum of three CpG sites and we only retained windows with among-sample variances greater than 0.05 to reduce computational burden when constructing the network (N ¼ 59,803 windows). We calculated connectivity between all pairs of regions using the bi-weight midcorrelation raised to the power of 18 to approximate a scale-free network as implemented in the WGCNA package in R (Langfelder and Horvath 2008). Modules of correlated methylation signatures were inferred using hierarchical clustering of the topological overlap dissimilarity matrix and a dynamic tree-cutting algorithm. The modules were constructed using a block-wise approach with a maximum of 30,000 regions allowed in each block and all blocks were then merged to form the final modules. The association of methylation modules with phenotypes was assessed by correlating module eigenvectors scores (the first axis of a principal component analysis conducted on all the regions within each module) with phenotypic values for each individual using the bi-weight midcorrelation. Module-traits correlations with P-values <0.05 were retained for further analysis. For each significant module-trait correlation we assessed whether the any regions within the module overlapped with previously identified DMRs found between SAS and wild fish. We assessed the statistical significance of these associations using a resampling procedure to compare the number of DMR overlaps with those based upon a random draw of regions for each module. We tested for differences in average module membership (region correlation with module eigenregion) and gene significance (region correlation with phenotype) between DMR-overlapping and non-DMRoverlapping regions using t-tests in R.
Wild salmon may experience selection during their time in the marine environment (i.e., smolt to adult; Bourret et al. 2014) where mortality can range from 65% to 99% (Chaput 2012). In contrast, SAS salmon may experience relaxed selection during this time in the hatchery, which likely explains the lower observed mortality rate for the studied cohort of 54% (from smolt collection to initiation of the juvenile experiment). To test for genetic divergence between SAS and wild fish we first used BisSNP v1.0.0 (Liu et al. 2012) to call single nucleotide polymorphisms from the aligned bisulfite sequencing reads. For this analysis, we retained all 16 adult samples and one individual per full-sib family from the juvenile dataset (N ¼ 8). Juvenile individuals were chosen to maximize the number of successfully called genotypes. We required a minimum depth of coverage of 8X to call an individual genotype and we retained only SNPs with a successful genotyping rate of 80% (N ¼ 24 individuals). We further removed SNPs with minimum allele frequency of 5% (minimum of 2 alternate alleles). We used AMOVA implemented in pegas v0.13 (Paradis 2010) to test for genome-wide differentiation. We used BayeScan v2.1 (Foll and Gaggiotti 2008) with a liberal prior odds setting of 10 as well as OutFLANK v0.2 (Whitlock and Lotterhos 2015) to identify potential outliers between SAS and wild groups. For the OutFLANK analysis we first built the genome-wide null distribution of divergence based on a linkage pruned set of 12,221 SNPs (obtained with "-indep 50 5 1" using PLINK v0.19) and tested for outliers in the whole dataset based on this distribution. We also employed a polygenic framework to test for subtle correlated changes across many alleles using redundancy analysis (Forester et al. 2016) implemented in the vegan v2.5-6 R package (Oksanen et al. 2019).

Adults: differential methylation in sperm
To identify the potential for intergenerational transmission of environmentally induced variation in DNA methylation acquired during gamete maturation, we quantified cytosine DNA methylation in sperm for over 16.4 million sites in the cytosinephosphate-guanine (CpG) context with greater than 5X coverage, but less than 20X coverage in 75% (9/12) of samples. Mean coverage per individual was 8.3X (range 4.6X-19X; Supplementary  Table S1). CpGs in salmon sperm were highly methylated (average methylation: 89%) and exhibited a bimodal distribution, with $5% of CpGs nearly devoid of methylation (<5%) and 91% of CpGs having methylation levels >80% (Supplementary Figure S1), consistent with levels reported in other salmonids (Gavery et al. 2018) and fish in general (Jiang et al. 2013;Potok et al. 2013;Ortega-Recalde et al. 2019;Skvortsova et al. 2019).
We identified differential methylation for individual CpGs between wild and SAS adults using beta-binomial models (Feng et al. 2014). There were 4998 differentially methylated cytosines (DMC; P-value < 0.001) identified between adult SAS and wild salmon sperm that were grouped into 284 differentially methylated regions (DMR; Figure 2A; Supplementary Table S2). Regions ranged in size from 51 to 2229 bp, contained between four and 34 CpGs, and comprised 90.2% of DMCs with FDR <0.05. The average magnitude of methylation difference between SAS and wild salmon for the identified DMRs was 39% (range: 8% to 70%). DMRs in SAS fish were twice as frequently hypo-methylated relative to wild fish than hyper-methylated and this difference was highly significant (68%: 193/284 hypo-DMRs; 32%: 91/284 hyper-DMRs; binomial P-value < 0.001). DMRs overlapped 237 genes or their cis-regulatory contexts (within 5000 bp of gene features; Supplementary Table S3). Gene ontology enrichment analysis revealed DMRs in sperm were associated with genes significantly enriched (P-value < 0.05) for a variety of functions in signal transduction pathways, brain development, tissue differentiation, muscle development and contraction, and chromatin silencing (Supplementary Table S4).
Juveniles: differential methylation in liver F1 SAS juveniles tended to be longer (SAS: 65.2 6 6.6 mm, wild: 63.2 6 6.6 mm; mean 6 SD; Figure 3A) and heavier (SAS: 3.42 6 1.0 g, wild: 3.14 6 1.0 g; mean 6 SD; Figure 3B) than F1 wild juveniles at 10-months of age, but these differences were dwarfed by among-family variation and were not statistically significant when controlling for the number of families investigated (length: . Two DMRs (C) were found in common between the two tissues and exhibited similar patterns of differential methylation in both adults and juveniles.
To investigate the potential for inherited methylation patterns in juvenile liver tissues, we quantified DNA methylation at over 23.1 million CpG sites covered with a minimum of 5X and a maximum of 20X in 75% (12/16) of F1 wild and SAS juveniles. Mean coverage per individual was 9.4X (range: 3.5X to 13.2X; Supplementary Table S1). In contrast to sperm, the more metabolically active liver tissues exhibited average methylation levels of approximately 80%. Juvenile liver tissue also exhibited a bimodal distribution of methylation where $5% of CpGs were unmethylated (<5%), 76% of CpGs with methylation >80%, and a larger fraction of sites with intermediate methylation levels (19% liver CpGs vs 4% sperm CpGs with methylation fraction between 5% and 80%; Supplementary Figure S2).
Differentially methylated CpGs (P-value < 0.001; N ¼ 5654) between wild and SAS juvenile offspring were organized into 346 DMRs that ranged in size from 51 to 2131 bp, contained between 4 to 40 CpGs, and covered 98% of DMCs with FDR < 0.05 ( Figure 2B; Supplementary Table S5). Similar to sperm cells, hypomethylation was almost twice as common in SAS juvenile liver tissues compared with hyper-methylation and this difference was highly significant (62%: 215/346 hypo-DMRs; 38%: 131/346 hyper-DMRs; binomial P-value < 0.001). The average magnitude of methylation differences between SAS and wild offspring were comparable (mean: 30%; range: 5%-52%). DMRs in juvenile liver tissues overlapped 274 genes or their cis-regulatory context (Supplementary Table S6). Over-represented biological functions of these genes reflected nervous system development and regulation, muscle development and contraction, signal transduction pathways, and immune system processes (Supplementary Table  S7).
Despite the fact that the adults we sequenced are not the parents of the juveniles we sequenced we found overlap of DMRs between the two tissues and life stages (2/622 total DMRs; Figure 2C) that was more than would have been expected by chance (1000 permutations; P < 0.001). The shared regions exhibited the same direction of differential methylation in both tissues (Supplementary Figure S3, A and B) and hierarchical clustering of these regions by individual largely recapitulated the SAS vs wild groupings ( Figure 2C). These regions are in proximity (<20 kb) to genes involved in immune response, tissue differentiation and organ development, and G protein-coupled receptor signaling (Supplementary Figure S3, A and B and Table S8). Five additional genes were overlapped by DMRs in both tissues but the DMRs in each tissue were not located at the same sites (Supplementary  Table S8). Of these, metabotropic glutamate receptor 4-like (GRM4-l) had regions that were in close proximity (<5000 bp). The DMRs in proximity to GRM4-l were hypo-methylated in adult SAS sperm and hyper-methylated DMR in juvenile SAS liver tissue (Supplementary Figure S3C).

Methylation influence on juvenile phenotype: comethylation network analysis
As a means to investigate whether methylation influences juvenile phenotypic variation (i.e., length, weight, and condition factor), we used a network-based approach (Langfelder and Horvath 2008) to identify modules of correlated methylation signatures across the juvenile liver tissue samples and tested for module associations with juvenile phenotypes. To construct the network, we first binned methylation in nonoverlapping 100 bp windows 63.2 6 6.6 mm) and heavier (B; SAS: 3.42 6 1.0 g, Wild: 3.14 6 1.0 g) than wild juveniles when reared in a common hatchery environment, wild juveniles had similar condition (C; SAS: 1.20 6 0.06, Wild: 1.21 6 0.06). All measurements are mean 6 standard deviation. Controlling for amongfamily variance using a linear model with a random-effect for family rendered none of the comparisons statistically significant (length: and selected the windows (i.e., regions) with among-individual variances greater than 0.05 (N ¼ 59,803 regions). We identified 124 modules that included a total of 4179 genomic regions. Eighteen modules exhibited significant correlations with at least one phenotype (P < 0.05; Supplementary Figure S4). Modules were enriched for a variety of signaling pathways and developmental processes relevant to the correlated phenotypes (e.g., growth factor signaling, skeletal muscle development; Supplementary Table  S6). The strongest association existed between the "navajowhite1" module and condition factor (r ¼ 0.69, P ¼ 0.003). The regions in this module were in proximity (5 kb) to genes enriched for mesoderm development and regulation of transcription (Supplementary Tables S5 and S6). In particular, the gene encoding insulin-like growth factor I (IGF-1), a hormone produced in the liver and secreted into the circulatory system, that is a key regulator of growth in muscle and skeletal tissues (Ohlsson et al. 2009) was associated with this module. Elevated methylation of IGF-1 in juvenile livers was associated with reduced weight and length, but better condition factor in juvenile salmon (Supplementary Figure S4). Three modules exhibited significant overlap with DMRs identified between juvenile SAS and wild salmon (permutation tests; P < 0.001). DMR-overlapping regions in these three modules had consistently higher module centrality than non-DMR-overlapping regions (Figure 4; purple: t 35.0 ¼ 2.9, P ¼ 0.007; yellow4: t 26.0 ¼ 3.0, P ¼ 0.006; darkorange: t 24.4 ¼ 2.6, P ¼ 0.01).
To account for the possibility of selection causing genetic differences between juvenile SAS and wild salmon and explaining the observed phenotypic differences, we also genotyped individuals for 974,219 single nucleotide polymorphisms (excluding CpG context C/T and A/G SNPs to avoid confounding methylation variation with allelic variation) using the aligned bisulfite sequencing reads and performed outlier tests. Overall, we failed to find support for a genome-wide average F ST larger than zero (AMOVA: 1000 permutations, P ¼ 0.77). Using outlier detection algorithms, we did not detect significant shifts in allele frequencies between SAS and wild salmon using BayeScan or a polygenic framework (RDA; R 2 ¼ 0, P ¼ 0.71), and we only identified two outlier SNPs with FDR < 0.01 in OutFLANK (Supplementary Figure S5 and Table S8), suggesting an overall absence of differential selection at the genome level between SAS and wild salmon.

Maturation environment influences gamete methylation
Our results demonstrate that environmental variation (i.e., growth and maturation in a natural vs hatchery setting) Enrichment of SAS-wild differentially methylated regions (DMR) in juvenile phenotype-associated methylation module cores (circles: red ¼ DMR overlap, gray ¼ no overlap). Three methylation modules correlated with juvenile phenotypes (A andD: purple-weight, r ¼ À0.34, P ¼ 0.04; B andE: yellow4-condition factor, r ¼ À0.23, P ¼ 0.04; C andF: darkorange-condition factor, r ¼ À0.53, P ¼ 0.03). All modules exhibited a significant correlation between module membership (x-axes; MM ¼ absolute value of the correlation of region methylation with the main axis of module variation) and gene significance (y-axes; GS ¼ absolute value of the correlation of region methylation with the phenotype). DMR overlapping regions in these modules were more highly correlated with the main axis of variation of the modules than non-DMR overlapping regions (A-purple: t 35.0 ¼ 2.9, P ¼ 0.007; B-yellow4: t 26.0 ¼ 3.0, P ¼ 0.006; C-darkorange: t 24.4 ¼ 2.6, P ¼ 0.01). DMR overlapping regions were also more centrally located and highly connected in module networks than non-DMR overlapping regions (D-F).
experienced after approximately 2 years of common rearing in their natural riverine environment alters both DNA methylation in salmon sperm as well as DNA methylation of hatcheryproduced progeny with SAS parents. The maturationenvironment effect we demonstrated strongly suggests, as other have reported (e.g., Dias and Ressler 2014;Sciamanna et al. 2019), that the Weismann Barrier is permeable and that information perceived by the soma can be incorporated into the germline via epigenetic mechanisms. Our results are consistent with recently published work indicating that DNA methylation patterns are dynamic in teleost gonads and sensitive to environmental variation (Todd et al. 2019). Several lines of evidence support a mechanism of environmentally mediated DNA methylation remodeling during gamete maturation. First, environmental differences in early life are known to influence gamete methylation (Gavery et al. 2018;Rodriguez Barreto et al. 2019). Second, multiple copies of DNA methyltransferase 3 (DMNT3), the methyltransferase responsible for the addition of new DNA methylation, have been retained following successive genome-duplication events (i.e., ohnologs) for salmonid and other teleost fishes (Liu et al. 2020). In Rainbow Trout (Oncorhynchus mykiss), certain DMNT3 ohnologs are expressed in spermatozoa during late spermatogenesis (i.e., a few weeks before spawning; Liu et al. 2020), thus providing a mechanism by which salmonids may alter DNA methylation in their gametes until days or weeks before spawning. Third, teleost fishes do not appear to experience genome-wide reprogramming of paternal methylation patterns following fertilization (Jiang et al. 2013;Potok et al. 2013) (Wang and Bhandari 2019), the presence of shared differentially methylated regions between our adult and juvenile datasets despite the fact they are not related suggests that at least some genomic regions escape reprograming even if it does also occur in salmon. Thus, adult salmon may be able to transmit heritable information to their offspring about the physical or biological environments they experience immediately prior to spawning. As such, our results are consistent with the hypothesis that epigenetic mechanisms can facilitate transgenerational plasticity (Bell and Hellmann 2019).
Transgenerational plasticity is theorized to evolve when environmental variability is sufficiently stable or predictable such that adults can transmit relevant information about the environment to their offspring (McNamara et al. 2016). If only early-life epigenetic signals were capable of being transmitted intergenerationally, transgenerational plasticity may not have been expected to evolve for salmon whose lives generally span both temporally and spatially diverse environments (e.g., freshwater habitats during early stages and marine habitats during late juvenile to adult stages), and whose life histories involve limited parental care (Thorstad et al. 2010). Our demonstration of the potential for adults to transmit environmental information acquired later in life to their offspring suggests transgenerational plasticity in salmon may be an important factor contributing to life-history variation and adaptive responses to environmental change.

Origins of hatchery-induced differential methylation
There is an unresolved question of whether epigenetic differences arising as a result of hatchery exposure occur due to deterministic processes (i.e., adaptive responses potentially arising from existing molecular machineries as a result of past selection) or stochastic processes (i.e., random environmental perturbations of wild-type methylation patterns). Several lines of evidence suggest the patterns we observed originate at least partly from deterministic processes. First, if methylation changes were completely stochastic, we would expect no bias in the direction of methylation changes. In contrast, our results show that differential methylation in SAS fish is strongly biased (about twofold) toward hypomethylation in both sperm and juvenile livers. Second, despite the fact that our juveniles and adults originated from different cohorts and that juvenile livers will have undergone tissuespecific methylation reprogramming from the state observed in sperm, we detected more DMRs in common between the datasets than expected by chance. Finally, the regions in phenotypically correlated methylation networks that overlap SAS vs wild DMRs are significantly biased toward being centrally located regions in the co-methylated networks providing further indications that these are not random associations.
More generally, prior research has established parallel signatures of DNA methylation divergence in response to early-life hatchery rearing in three populations of Coho Salmon (Oncorhynchus kisutch; Le Luyer et al. 2017) and hatchery-reared Rainbow Trout exhibit a significant proportion (approximately 20%) of hatchery-origin DMRs that are shared between red blood cells and sperm (Gavery et al. 2018). There are similarities (discussed in detail below) in the biological functions impacted by the epigenetic signatures of response to early-life hatchery rearing for Rainbow Trout (Gavery et al. 2018 (Konstantinidis et al. 2020). We note that the design of our study cannot directly attribute inherited changes solely to the paternal methylation pattern as observed in Zebrafish (e.g., Skvortsova et al. 2019). In particular, maternal effects driving variation in DNA methylation have been described (Venney et al. 2021) and thus future work should consider the dual roles of maternal and paternal origin in methylation inheritance of hatchery reared salmon.

Epigenetic signatures of domestication
We identified seven genes associated with SAS versus wild differential methylation that have been reported in previous studies of hatchery induced differential methylation in salmonids. Phosphatidylinositol 3-kinase regulatory subunit alpha (PIK3R1) was differentially methylated in SAS sperm cells and has previously been identified as differentially methylated in sperm of Atlantic Salmon reared in a hatchery from fertilization (Rodriguez Barreto et al. 2019). This gene is an important regulator at the center of many growth factor and hormone signaling pathways (Cantley 2002) and has also previously been identified as potentially target of selection in domesticated salmon (Liu et al. 2017). Three additional genes reported in previous studies that we report as differentially methylated in sperm cells all have roles in the nervous system including growth and differentiation (NRG2; Britsch 2007), cell-cell adhesion (PCDHGC5; Wang et al. 2002), and neurotransmitter release (STXBP5L; Geerts et al. 2015).
The remaining three genes we detected in common with other studies were all differentially methylated in juvenile liver tissue with reported functions in other organisms of epidermal tissue development (BCR; Dubash et al. 2013), cell adhesion and cytoskeleton organization during neuron differentiation (CTNNA2; Schaffer et al. 2018), and neuron differentiation (ARHGAP32; Nakamura et al. 2002). BCR has previously been identified as differentially methylated in hatchery-origin Coho Salmon white muscle (Le Luyer et al. 2017), CTNNA2 in both red blood cells and sperm of hatchery-origin Rainbow Trout (Gavery et al. 2018), and ARHGAP32 in sperm of Atlantic Salmon (Rodriguez Barreto et al. 2019). In general, patterns of DNA methylation across studies implicate regulation of cell differentiation and developmental processes with a particular enrichment of genes involved in neuron differentiation.
We identified several glutamate receptors as being differentially methylated either in sperm (GRIK5 and GRM4) or liver (GRM4, GRID2, and GRM3). Glutamate receptors are well-known targets of selection in many domestic animals (O'Rourke and Boeckx 2020) and in domestic Atlantic Salmon specifically, where apparent positive selection has been reported for structural variants (i.e., insertions, deletions, and inversions) associated with synaptic genes, including multiple glutamate receptors (Bertolotti et al. 2020). Glutamate receptors are also differentially methylated between wild and domestic European Sea Bass (Anastasiadi and Piferrer 2019). Glutamatergic signaling is an important excitory driver of the hypothalamic-pituitary-adrenal (HPA) axis that, among other functions, mediates organismal stress responses and aggression. It has recently been hypothesized that selection on these pathways underlies "tameness" in domesticated animals and attenuated stress responses under crowded conditions (O'Rourke and Boeckx 2020). This raises the hypothesis that adult SAS fish could transmit information about the crowding they experienced prior to spawning to their offspring in order to prime them for a highly competitive environment upon hatching (Christie et al. 2012).
In summary, were the epigenetic effects induced by hatchery environments truly random, it would be very unlikely to detect particular genes or pathways across multiple studies as revealed by the available literature. Altogether, the compiled evidence supports the hypothesis that there is a certain degree of conservation in the DNA methylation changes in response to captive rearing across a broader taxonomy of teleost fishes.

Consequences of hatchery-induced methylation for offspring phenotypes
We identified several correlated methylation profiles that were statistically associated with offspring phenotypes. Conceptually, this analysis identifies pathways or biological functions that are co-regulated by methylation. Of the 18 methylation modules that were correlated with juvenile phenotypes, several comprised genomic regions occurring in proximity to genes enriched for signaling pathways (i.e., muscle growth and differentiation, skeletal development, neural development, and immune system processes) directly relevant to the phenotypes being studied. In particular, the methylation module in juvenile livers with the strongest phenotypic correlation (i.e., navajowhite1) contained regions overlapping IGF-1 which is a hormone produced and released from the liver in response to growth hormone (GH) signaling that plays a key endocrine role mediating growth and differentiation of muscle and skeletal tissues and therefore body size (Ohlsson et al. 2009). Body size and condition factor are traits closely linked with juvenile salmonid survival and fitness (Quinn and Peterson 1996;Einum and Fleming 2000) and the GH/IGF-1 axis has been implicated in acclimation to saltwater (McCormick 2001), which is a major selective barrier for anadromous salmonids and has been identified as a deficiency of hatchery-reared fish (Shrimpton et al. 1994). Furthermore, relaxed selection on IGF-1 has been reported for landlocked populations of Atlantic Salmon where fish spend their entire life in freshwater (Kjaerner-Semb et al. 2021), a situation not unlike our hatchery-rearing conditions. Whether the differential methylation associated with IGF-1, we observed also reflects relaxed selection from not experiencing the marine environment, or reflects a directed change of methylation state to reach a fitness optimum for life in freshwater, or a hatchery environment, remains unclear. More work is required to understand the relative contributions of different aspects of smolt to adult rearing conditions on variance in DNA methylation and then ultimately link them to more complex phenotypes including fitness. As a start, our results demonstrate that differences in the DNA methylation state of this important growth-regulating gene have the potential to exert influence on salmon growth and developmental trajectories that are likely to have real consequences for individual fitness.
Hatchery-induced differential methylation may directly influence both the fitness-related traits we quantified here and other more complex behavioral traits with fitness consequences at later stages of development than we have studied. Hatchery-induced DMRs overlapped regions centrally located in co-methylated modules that were associated with biological functions involved in brain neuron differentiation as well as the glutamate receptor GRM4 (i.e., module yellow4) suggesting that hatchery-environment-induced methylation-mediated behavioral changes (e.g., attenuated stress response to crowding; O'Rourke and Boeckx 2020) could have consequences for the growth trajectories of offspring. We also identified modules (i.e., purple) and differential methylation of genes not included in methylation modules that lie downstream of genes in important phenotype-associated modules (i.e., IGF-1 signaling pathway). Phosphatidylinositol-mediated signaling (i.e., module purple) and PI3KR1 in particular plays a key role in modulating the response to IGF-1 stimulation (Hakuno and Takahashi 2018) and thus hatchery-induced differential methylation may influence or modulate growth trajectories at a later step in that signaling cascade that, at least in part, may explain the significant differences we observed in phenotypes between SAS and wild juveniles. Collectively, our results suggest that hatchery-associated effects could indeed be mediated through DNA methylation with consequences for aspects of fish phenotypes with known relationships to their fitness.

Comparison of hatchery rearing approaches
Despite detecting DNA methylation differences between SAS and wild fish, our work also reveals some fundamental differences between the effects of early-rearing and later-life hatchery exposure for salmon. SAS fish (both SAS adults and their progeny) exhibited more hypo-methylation relative to wild fish (for both adult males and their progeny reared in a common environment), in contrast to previous work that has demonstrated predominantly hypermethylation of fish produced and reared in hatchery from the egg stage (Le Luyer et al. 2017;Gavery et al. 2018;Rodriguez Barreto et al. 2019). Reduced representation bisulfite sequencing (RRBS) in Coho Salmon and Rainbow Trout reported differential methylation at between 0.03% and 0.1% of analyzed sites or regions. Using similar criteria, our results showed differential methylation affected an order of magnitude less CpG sites (0.004%). This suggests that the potential for DNA methylation-mediated domestication effects caused by later-life hatchery exposure may not be as severe as those observed for salmon that experience early-life hatchery rearing. RRBS approaches are believed to preferentially target gene regulatory relevant regions of the genome (e.g., Anastasiadi and Piferrer 2019) and thus, because of a difference in techniques between studies (RRBS vs whole-genome bisulfite sequencing), our comparisons may be biased. However, the genomic distribution and density of regulatory relevant CpGs in nonmammalian vertebrates fundamentally differs from that of mammals (Long et al. 2013). Bioinformatic interrogation of our data indicates RRBS applied to our study would have assayed approximately 7 million CpGs and only detected 10% of the observed DMCs implying the above comparison is reasonably unbiased. Furthermore, it suggests that in salmon, and possibly fishes more generally, RRBS approaches fail to capture a significant proportion of biologically relevant methylation differences.
Like previous studies, we have demonstrated a lack of genome-wide sequence differentiation between hatchery-reared and wild fish (Christie et al. 2016;Le Luyer et al. 2017;Gavery et al. 2018). This result suggests that the selective regime imposed by the hatchery environment over one generation was not strong enough to cause widespread differentiation. In turn, it also suggests that despite the high levels of mortality during the marine phase of salmons' lives (65-99%; Chaput 2012), selection in the marine environment may not be important enough to cause widespread, temporally consistent, changes in allele frequencies between wild and SAS salmon. Previous work in Atlantic Salmon has reported consistent allele frequency changes over the marine migration period for only one of two populations studied indicating patterns of differentiation due to the marine environment are spatially and temporally variable (Bourret et al. 2014). It is difficult to know if the two outliers we detected result from selection in the hatchery or marine environments. In spite of this, our results clearly implicate a stronger role for epigenetic factors and not differences in genetic variation in hatchery-related phenotypic divergence.
We have demonstrated the potential for environmental effects to be propagated to offspring for salmon who experience hatchery environments during maturation via the intergeneration transmission of DNA methylation. Our experiments on juvenile fish were conducted in a laboratory setting and thus whether these effects are also detectable and have fitness consequences in the true context of SAS program where SAS individuals are released and reproduce naturally in the wild remains unknown. Genotype-by-environment interactions are pervasive in salmonids (Vandersteen et al. 2019) and so are epigenotype-byenvironment interactions (Christensen et al. 2021). As such, there is an urgent need to evaluate the interaction between SAS and wild rearing on offspring DNA methylation and development in a natural environment. Other sources of epigenetic information (i.e., small RNAs) are also well known to mediate intergenerational effects (Sciamanna et al. 2019) that may well mediate the phenotypic differences between individuals we have observed. On the other hand, multiple epigenetic mechanisms often function together to affect phenotypic changes (Cavalli and Heard 2019) and future work unraveling the mechanistic basis of hatchery-induced phenotypic effects will need to clarify the potential contributions of other epigenetic mechanisms, their relative importance, and the degree to which these effects are reversible following the cessation of the environmental exposure. Only then will the evolutionary consequences of environmentally induced epigenetic variation in these systems be globally understood.