Altered spawning seasons of Atlantic salmon broodstock genetically 1 and epigenetically influence cell cycle and lipid-mediated regulations 2 in their offspring 3

17 Manipulating spawning seasons of Atlantic salmon ( Salmo salar ) is a common practice to facilitate 18 year-round harvesting in salmon aquaculture. This process involves adjusting water temperature and 19 light regime to control female broodstock maturation. However, recent studies have indicated that 20 altered spawning seasons can significantly affect the nutritional status and growth performance of 21 the offspring. Therefore, gaining a deeper understanding of the biological regulations influenced by 22 these alterations is crucial to enhance the growth performance of fish over multiple generations. 23


Introduction
The Atlantic salmon (Salmo salar) follows an anadromous life cycle, beginning in freshwater, migrating to seawater for most of its life, and returning to freshwater to spawn.In salmon aquaculture, the timing of spawning is adjusted to ensure a year-round availability for harvesting.This manipulation involves controlling abiotic factors such as water temperature and photoperiod to influence the spawning timing [1,2].Additionally, apart from the traditional open sea-based cages, a land-based approach utilising recirculating aquaculture systems (RAS) is employed to replicate the salmon's natural lifecycle more accurately, encompassing both saltwater gonad maturation and final freshwater maturation [3].Despite the ongoing and successful practices of spawning season alterations in salmon aquaculture, recent studies have reported that these practices can impact the nutritional status of broodstock, leading to negative effects on the growth performance in the next generation [4,5].
Nutritional programming, in the context of fish nutrition, refers to strategically manipulating the feeding regime of fish during critical developmental stages to enhance their long-term growth performance [6][7][8].It can be expanded into intergenerational programming by providing specific nutrients or dietary components to the parents before reproduction [8,9].Hence, spawning season alterations of salmon broodstock can be seen as an example of intergenerational nutritional programming even though the nutritional manipulation occurs indirectly through abiotic factors.The principal target of the programming is to ensure the production of healthy and robust offspring in this context.
While various nutrients influence the health condition and growth performance of fish, including Atlantic salmon [10][11][12], altered spawning seasons have appeared to affect the status of multiple metabolites and nutrients including micronutrients, amino acids, and lipid classes [4,5].These nutrients are also closely associated with vital biological functions essential for growth, such as one carbon (1C) metabolism, the citric acid cycle, and the Cahill cycle.Key micronutrients like vitamin B12, folate, vitamin B6, and methionine are integrated with the 1C metabolism [13,14], which also regulates cellular methylation potential through S-adenosylmethionine (SAM) and S-adenosylhomocysteine (SAH), critical for various cellular processes [15].The citric acid cycle enables the catabolism of carbohydrates, fats, and amino acids, releasing stored energy crucial for growth and development [16], while the Cahill cycle facilitates the transport of amino groups and carbons between muscle and liver [17], both playing essential roles in various developmental processes.Amino acids are essential for growth and tissue development, in addition to being the building blocks of proteins [18], while lipids serve as major energy reserves and contribute to the formation of cell membranes, also supporting various developmental processes [19].
Beyond the importance of nutrients, understanding genetic influences within the context of nutritional programming provides invaluable insights from a molecular biology perspective.Recent studies have revealed that alterations in micronutrient levels within the feed can have genetic effects on various biological pathways in salmon, including lipid metabolism, protein synthesis, and post-transcriptional regulation [20,21].Moreover, an additional layer of information comes from the analysis of DNA methylation, which offers valuable epigenetic perspectives.DNA methylation, being an epigenetic regulatory process, responds to environmental stimuli, governing cell stability, differentiation, and development [22,23].As DNA methylation is inherited during cell replication, it can lead to changes in epigenetic regulation over the long term [22,24].
The present study investigates the genetic and epigenetic effects resulting from four different spawning seasons.These seasons include the natural sea-pen-based spawning in November, an accelerated RAS-based spawning occurring five months earlier, and two variations of sea-pen-based spawning: one with a two-month advance and another with a two-month delay relative to the natural spawning season.The main analysis of the study focuses on liver samples collected from start-feeding larvae at the end of the endogenous feeding phase.The liver serves as a key target organ for various biological reactions crucial to growth, and using start-feeding larvae helps eliminate the influences of exogenous feeding.
To uncover both genetic and epigenetic impacts comprehensively, we employed RNAsequencing (RNA-seq) for gene expression and reduced-representation bisulfite sequencing (RRBS) for DNA methylation analysis.Our findings provide valuable insights into the significance of the nutritional state in broodstock, as it can affect the growth of their offspring.Additionally, our study offers potential usages of gene expression and DNA methylation data as powerful tools for assessing current and future growth performance, which can significantly benefit a wide range of nutritional programming studies.

Ethical considerations
The experiment adhered to the ARRIVE guidelines for design and reporting.Broodstock and embryos used in the study were obtained from AquaGen's commercial production facility of Atlantic salmon at their breeding station in Kyrksaeterøra, Norway.
The conditions and protocols applied in this experiment were identical to those used in commercial production, including certain confidential information protected under commercial interests.Throughout the sampling procedures for broodstock and larvae, anesthetics were administered by the supplier's instructions and compliance with both Norwegian and European legislation on animal research.
Formal approval from the Norwegian Animal Research Authority (NARA) was not required for this experiment.This exemption is because the experimental conditions fall under practices recognized as commercial animal husbandry, thereby exempting them from the European Convention on the protection of animals used for scientific purposes (2010/63/EU), specifically under Article 5d.Additionally, the Norwegian ethics board also approved the experiment by the Norwegian regulation on animal experimentation, § 2, 5a, d, for activities classified as "non-experimental husbandry (agriculture or aquaculture)" and "procedures in normal/common breeding and husbandry".

Experimental design and sampling
As the details of experimental design and sampling methods have been described elsewhere [4,5], we only present summarized descriptions here.The study utilized broodstock and embryos obtained from AquaGen's breeding station in Kyrksaeterøra, Norway.The female broodstock, sourced from the 15th generation, covered four spawning seasons: off-season, early, normal, and late seasons.
The off-season broodstock were reared in a land-based recirculating aquaculture system (RAS) in saltwater (12% salinity), while the rest of the broodstock were raised in seabased net pens.All broodstock were fed to satiation with broodstock feed (EWOS Breed 3500) until transferred to freshwater.The normal season broodstock matured naturally under ambient photoperiod and temperatures, whereas artificial controls were applied to regulate maturation for the other spawning seasons.Specifically, the off-season and early season broodstock were subjected to five-and two-month advance maturation, respectively, while the late season broodstock were subjected to a two-month delay.The fish received no further feeding after the transfer to freshwater.The maturation process in freshwater lasted for 109, 163, 120, and 166 days for off-season, early, normal, and late season broodstock, respectively, before fertilization to obtain the next generation (Supplementary Table S1).
From each spawning group, five females from separate tanks were sacrificed using Benzoak (200 mg/L, ACD Pharmaceuticals AS).After stripping, each female's growth measures, including body weights (n=5 per spawning season), were recorded.Liver samples were dissected and flash frozen in liquid N2 for nutrient analysis.All oocytes were fertilised with cryopreserved sperm from two pooled males.Newly fertilized eggs (2 d°) were flashfrozen in liquid N2 for nutrient analysis (Supplementary Table S1).
At the end of the endogenous feeding phase (start-feeding larvae 979-994 d°; Supplementary Table S1), 36 larvae from each broodstock were randomly selected, euthanised with an overdose of buffered tricane methanesulfonate (MS222, Pharmaq, Norway), and their growth measures, including body weights (n=180 per spawning season), were recorded.Additionally, liver samples were dissected from larvae and flash frozen in liquid N2 for RNA and DNA analysis.

Growth measure and nutrient analysis
We collected the data of body weights for broodstock and larvae from previous studies [4,5].
In addition, we obtained the data on selected nutrients and metabolites from the same studies, specifically vitamin B6 & B12, folate, S-adenosylmethionine (SAM), Sadenosylhomocysteine (SAH), cholesterol, the sum of six different lipids, glutamate, Lserine, L-lysine, L-alanine, L-glutamine, urea, and B-alanine [4,5].To analyse the data from all four spawning seasons, we performed ANOVA followed by Tukey's HSD test.

DNA and RNA extraction
RNA extraction was carried out using the BioRobot EZ1 and EZ1 RNA Universal Tissue kit (Qiagen), followed by DNase treatment using the Ambion DNA-free DNA removal kit (Invitrogen, USA) according to their respective protocols.The RNA quantity was assessed using the NanoDrop ND-1000 Spectrophotometer (Nanodrop Technologies) and the Agilent 2100 Bioanalyzer with the RNA 6000 Nano LabChip kit (Agilent Technologies).For DNA isolation, the DNeasy Blood & Tissue Kit (Qiagen, Cat.No. #69506) was used as per the manufacturer's protocol.Liver samples underwent RNase A treatment (provided by the Qiagen kit, 50ng/µL, 10 min at room temperature) followed by proteinase K treatment (New England Biolabs, #8102S 20µg/µL, 1.5 h at 55°C).DNA was eluted in Milli Q water.DNA quantification was performed using the Qubit High Sensitivity Assay (Life Technologies #Q32854).Detailed methods for both DNA and RNA extraction were described elsewhere [20].

Library preparation of high throughput sequencing
Liver samples were sequenced at two different sequencing facilities for specific analyses.For RNA-seq, the liver samples were processed at the DeepSeq sequencing facility at Nord University, Bodø, Norway, and the library preparation was conducted using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs).Sequencing was performed on the NextSeq500 machine (Illumina).For RRBS, the liver samples were processed at the CeMM Biomedical Sequencing Facility, Vienna, Austria.The genomic DNA was extracted, digested by MspI, and bisulfite-converted before library preparation.The RRBS libraries were sequenced on Illumina HiSeq 3000/4000 instruments.Detailed methods of the library preparation for RNA-seq and RRBS were described elsewhere [20].

Atlantic salmon genome and genomic annotation
The reference genome (ICSASG version 2) and RefSeq data (version 100) for gene annotation were downloaded from the NCBI web site (https://www.ncbi.nlm.nih.gov/assembly/GCF_000233375.1).In cases where gene symbols were either outdated or unavailable from RefSeq data (version 100), gene symbols from a newer version of RefSeq data (version 102) and UniProt [25] were used.The genome regions were annotated into intron, exon, three promoter regions, and flanking regions.Promoter regions were divided into three categories based on their distance from transcription start sites: P250 (1-250 bp), P1K (251-1000 bp), and P5K (1001-5000 bp).Flanking regions were defined as 5000 upstream from P5K and 10000 downstream from the 5' end of the gene.

Pre-processing of high throughput sequencing
For both RNA-seq and RRBS data, the same pre-processing procedures were followed as described in a previous study [20].This involved trimming the reads using Cutadapt [26] and Trim Galore! (Barbraham Institute), aligning the trimmed reads to the reference genome using STAR [27] for RNA-seq and Bismark/Bowtie 1 [28,29] for RRBS with their default parameters.The mapped RNA-seq reads were quantified using featureCounts [30], while the mapped RRBS reads were processed by Bismark [28] for methylation calling and CpG site extraction.The samples were divided into four spawning season groups for further analysis: off-season, early, normal, and late seasons, and clustering analysis was performed using the factoextra package (https://CRAN.R-project.org/package=factoextra).RNA-seq counts were subjected to a variance stabilizing transformation (VST) using DESeq2 [31] before principal component analysis (PCA).
Raw reading data were uploaded to the repository of the Sequence Read Archive (SRA) on the NCBI web site and available under the accession numbers PRJNA680425 and PRJNA642998 for RNA-seq and RRBS, respectively.

Differential gene expression analysis
Differentially expressed genes (DEGs) were identified using DESeq2 [31] with an adjusted pvalue cut-off of less than 0.05.P-values were adjusted by the Benjamini-Hochberg procedure [32].From the candidate DEGs, those with absolute log fold changes (LFCs) greater than or equal to 1.2 were selected as the final set of DEGs by using the lfcThreshold argument of the results function of DESeq2 [31].Shrunken LFCs values, calculated by the normal shrinkage method provided by DESeq2 [31], were used for both heatmaps showing functional annotation results and scatter plots showing merged results between DNA methylation and gene expression differences.

Functional annotation with KEGG
Over-representation analysis (ORA) was performed on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [33] using clusterProfiler [34].Gene lists of DEGs identified from three pair-wise comparisons were used as input, and enriched pathways were determined based on adjusted p-values less than 0.05.P-values were adjusted by the Benjamini-Hochberg procedure [32].

Methylation rate analysis
Mapped RRBS reads were filtered by methylKit [35], discarding reads with coverage less than or equal to 10 and above the 99.9th percentile.Additional RRBS data from two other studies were obtained for post-smolt and harvesting stages [20,21] and underwent the same pre-processing procedures to compare the distributions of methylation rates.Raw data of these RRBS sequencing data can be found in SRA on the NCBI site under the accession numbers of PRJNA680423 and PRJNA628740 for post-smolt and harvesting stages, respectively.Ridge plots were generated using the ggridges R library (https://CRAN.Rproject.org/package=ggridges).

Differential methylation analysis
Differentially methylated CpG sites (DMCs) were identified by methylKit [35], and Q-values were calculated using the sliding linear model (SLIM) method [36].CpGs were identified as DMCs when Q-values were less than 0.01 and methylation differences were greater than or equal to 15%.
To filter out DMCs that were strongly affected by altered spawning seasons, we used the counts of DMCs per gene in two promoter regions: P250 and P1K.Genes were sorted by the count of DMCs in a descendant order, separately for P250 and P1K, and genes that contained at least three DMCs were selected for further analysis.

Bioinformatics analysis
In-house R and Python scripts with Snakemake [37] were utilized for high-throughput sequence analysis, basic statistical analysis, and figure generation.

Altered spawning seasons impacted the weight and nutritional status of offspring
The present study examines the effects on gene expression and DNA methylation patterns in the liver of offspring hatched from four different spawning seasons: an off-season based on a recirculating aquaculture system (RAS) and three sea-pen based seasons (early, normal, and late) (Fig. 1a).Off-season spawning occurred in June, early season in September, and late season in January, relative to the normal spawning months of November.(see detailed dates and degree days in Supplementary Table S1).As the growth performance and nutritional status of the broodstock and offspring have already been assessed and discussed elsewhere [4,5], this section only provides a summary emphasizing the weight and the status of selected nutrients and metabolites.Broodstock weights showed no significant differences at spawning (upper part of Fig. 1b), but offspring weights at the larvae stage (979 and 994 degree days) showed noticeable variation (lower part of Fig. 1b).Specifically, off-season and early season larvae had lower weights compared to normal and late seasons.This trend was consistent with the observations in egg sizes, where eggs from off-season and early-season exhibited smaller sizes compared to the other two seasons (Supplementary Table S2).
Altered spawning seasons widely impacted the nutritional status of offspring eggs, particularly regarding 1C metabolism, lipid classes, citric acid cycle, and Cahill cycle [4,5], but with inconsistent patterns across nutrients and metabolites (summarised in Fig. 1c and Supplementary Tables S3-S6).Previous studies also reported that altered spawning seasons affected the nutritional status of the broodstock [4,5].Thus, altered spawning seasons significantly influenced the nutritional status of broodstock without impacting their weights, subsequently affecting offspring nutritional status and growth potential of their offspring.
To investigate the underlying genomic and epigenetic regulations associated with the observed differences in nutritional status and growth performance, we utilised RNA-seq for gene expression analysis and RRBS for DNA methylation analysis on liver samples collected from larvae at the end of the endogenous feeding phase.

Early and late spawning seasons exhibited similar gene expression patterns when compared to the normal season
We performed gene expression analysis on 20 liver samples divided into four groups (n=5) to investigate the impact of different spawning seasons to the offspring when compared to the normal season (detailed statistics for each sample are provided in Supplementary Table S7).
Principal component analysis (PCA) differentiated the off-season from the normal season when considering all mapped genes (left part of Fig. 2a).However, when focusing on the top 500 genes with high variances, the early season moderately deviated from the normal season (right part of Fig. 2a).In both cases, the late season fell between the normal and early seasons.In three pairwise comparisons for differential expression analysis (off-season vs. normal, early vs. normal, and late vs. normal), we identified 173, 337, and 265 differentially expressed genes (DEGs), respectively (Fig. 2b).All comparisons revealed a higher number of down-regulated genes than up-regulated genes.Additionally, Venn diagrams indicated that the DEGs from the off-season vs. normal comparison had very few overlapping genes (<1%) with the other two comparisons (depicted in the first and second diagrams in Fig. 2c).On the other hand, early vs. normal and late vs. normal comparisons shared approximately 25% of the DEGs (depicted in the third diagram in Fig. 2c).
The results from PCA and differential expression analysis suggest that the expression pattern of the off-season differed from both the early and late seasons when compared to the normal season.Furthermore, the early and late seasons exhibited similar gene expression patterns, with the late seasons showing more resemblance to the normal seasons.

Altered spawning seasons affected various biological pathways, including those related to metabolism, cellular processes, and organismal systems
To explore the impact of altered spawning seasons on biological pathways, we conducted functional analysis on DEGs using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.Over-representation analysis (ORA) identified a total of 11 enriched KEGG pathways, which were classified into three functional categories: metabolism, cellular processes, and organismal systems (Fig. 3).In the off-season vs. normal comparison, ORA identified only one enriched pathway, carbon metabolism (Supplementary Table S8).All five DEGs associated with carbon metabolism showed down-regulation (Fig. 3).
For the early vs. normal comparison, ORA revealed six enriched pathways: one in metabolism, three in cellular processes, and two in organismal systems (Supplementary Table S9).Metabolism-related pathways exhibited both up-and down-regulation, while most genes in cellular processes showed significant down-regulation.In organismal systems, PPAR signalling showed up-regulation, while progesterone-mediated oocyte maturation showed down-regulation (Fig. 3).
The late vs. normal comparison identified nine enriched pathways: three in metabolism, four in cellular processes, and two in organismal systems (Supplementary Table S10).Remarkably, this comparison shared five enriched pathways with the early vs. normal comparison, indicating similar impacts on biological pathways between them.Moreover, these pathways displayed consistent expression patterns in terms of up-and down-regulation between early and late seasons.Notably, genes associated with cellular processes exhibited strong down-regulation for both seasons (Fig. 3).
Altered spawning seasons affected multiple pathways in metabolism, cellular processes, and organismal systems.The off-season appeared to down-regulate carbon metabolism, while the early and late seasons primarily had strong suppressive effects on pathways associated with cellular processes.These findings may be linked to the reduced weight of off-season and early season larvae, as well as variations in the nutritional status of the eggs.

Altered spawning seasons had an impact on the expression of genes in a similar manner
To further investigate the genes commonly affected by multiple spawning seasons, we used overlapping DEGs (represented as intersections of Venn diagrams in Fig. 2c) by applying two filtering steps: selecting DEGs found in at least two comparisons and choosing the top five DEGs from each comparison based on adjusted p-values (Supplementary Table S11).This filtering approach revealed a total of 13 genes (Table 1).Among them, two genes (gene IDs: 106612264 and 106601246) lacked corresponding gene names and symbols.All genes, except pc and pitpnm2, consistently displayed up-or down-regulation across all spawning seasons, suggesting that all altered spawning seasons had a similar effect on this specific set of genes (Table 1).These genes can be classified into three categories by relevant biological functions: metabolism (pcyt1b, pc, klf15, and cyp7a1) [38][39][40], tissue development (prox3, crim1, and krt8) [41][42][43][44], and information processing, which can be further classified into three sub categories as signalling (bmp5 and pitpnm2) [38,45], transcription (tcf12) [46], and posttranslation (rhobtb4) [38].As these categories are typically associated with growth, all three altered spawning seasons affected a certain set of genes associated with growth regulation, despite PCA plots (Fig. 2a) and Venn diagrams (Fig. 2c) showing a distinct gene expression pattern of the off-season.

Overall methylation rates increased until the harvesting stage, except around transcription start sites
We conducted DNA methylation analysis on 20 liver samples from larvae across four different spawning seasons (n=5; detailed statistics in Supplementary Table S12) to investigate the epigenetic effects of altered spawning seasons on DNA methylation.To examine regional DNA methylation patterns, CpG sites were categorised into three genomic regions based on their positions relative to closely located genes: gene body (GB; including exons and introns), promoter (P; including P250, P1K, and P5K), and flanking regions (Fig. 4a; see Methods and Materials for definitions).Before analysing the effects of altered spawning seasons, we compared RRBS data from two other studies focusing on liver samples from different developmental stages: postsmolt (~31 weeks) and harvest (~45 weeks) [20,21].Distributions of methylation rates showed significant differences among life stages, except for P250 (Kolmogorov-Smirnov test; Supplementary Table S13).Ridge density plots displayed peak values ranging from 75% to 100%, excluding P250 and P1K (Fig. 4b).Notably, the peaks of the harvesting stage shifted closer to 100% compared to the larvae and post-smolt stages, supporting the trend of average methylation rates: larvae < post-smolt < harvest (average methylation rates are shown in text labels of Fig. 4b).This comparison suggests that methylation rates generally

Differential methylation patterns indicated higher dynamic regulation in promoter regions compared to other regions
We conducted six pairwise comparisons to identify CpG sites with differential methylation (Fig. 5a).Our objective was to capture a comprehensive pattern of methylation differences by considering all possible pair-wise combinations.Remarkably, all comparisons consistently demonstrated similar characteristics, for both the number of differentially methylated CpGs (DMCs) and the balanced ratio of hypo-and hyper-methylation (Fig. 5a).Further investigation focused on a set of 10,308 common DMCs, identified in at least one comparison.Distributions of methylation rates for these common DMCs exhibited a distinct left-skewed pattern across the exon, intron, P5K, and flanking regions (Fig. 5b).This observation strongly suggested that methylation differences resulting from altered spawning seasons primarily occurred within regions that were already highly methylated.In contrast, the P250 and P1K regions displayed nearly uniform methylation rates across a range of 0-100%, indicating significant methylation differences occurred irrespective of the underlying methylation rates.These findings imply a potentially greater level of dynamic regulation of DNA methylation in P250 and P1K compared to the other regions.

Filtering multiple DMCs identified key genes affected by spawning season alterations
We examined genes with multiple DMCs in three pairwise comparisons against the normal season.Since most genes had only a single DMC, especially in promoters (~75% in P250 and ~70% in P1K, Supplementary Fig. S2), we avoided identifying differentially methylated regions (DMRs) using a sliding window approach.Instead, we filtered out genes strongly supported by multiple DMCs based on specific criteria (see Methods and Materials), resulting in a total of 11 genes (Table 2), including six genes for off-season vs. normal and three genes each for the other two comparisons (Supplementary Table S17).exhibited either hypo-methylated or hyper-methylated DMCs, but not both, within a single comparison (Table 2).Furthermore, approximately 60% of the genes (7 out of 11) consistently displayed either hypo-methylation or hyper-methylation even across multiple comparisons ('All hypo' and 'All hyper' in the pattern column of Table 2, respectively).
Literature analysis associated these genes with two biological functions: development (inf2) and information processing, which can be further classified into four sub categories as transporter (slc31a1 and scp-2) [47,48], signalling (stk32a and lrrfip2) [25,49], post-translation (klhl11 and ttll12) [25,50], and transcription (helz2 and prrx2) [51,52].However, ctl2b and an uncharacterised gene (Gene ID: 106573279) had no clear associations with specific biological functions.Although these genes were potentially linked to growth performance, none of the 11 genes were identified as DEGs in our gene expression analysis.

Altered spawning seasons influenced gene expression and DNA methylation of cell cycle-related genes
To investigate potential associations between DNA methylation and gene expression, we focused on DMCs identified in three pairwise comparisons against the normal season (offseason vs. normal, early vs. normal, and late vs. normal) in P250 and P1K.By correlating their methylation differences with corresponding log-fold changes (LFCs), we found a total of 11 DMCs: four DMCs in P250 and seven DMCs in P1K associated with DEGs (Fig. 5c, Supplementary Table S18).Surprisingly, CpG sites with methylation differences greater than 30% lacked clear association with differentially expressed genes (light grey dots labelled as non-DEGs in Fig. 5c), suggesting that differentially methylated sites showed inconsistent links with differential gene expression.
These 11 DMCs were associated with seven genes that were DEGs in the corresponding comparisons (Table 3).We expanded this list of seven genes with additional information from all pair-wise comparisons (additional entries were marked with parentheses in Table 3).Among them, three genes (caprin-1, kifc1, and aurkb) exhibited consistent regulations of DMCs and DEGs compared to the normal season.Specifically, caprin-1 and kifc1 had hypo-methylated CpG sites, while aurkb had hyper-methylated CpG sites; all three genes showed down-regulated gene expression (Table 3).Notably, all three genes (caprin-1, kifc1, and aurkb) were associated with cell cycle regulation [53][54][55], while the remaining genes were linked to metabolism (cyp8b1, adrenodoxin, and lpin1) [33,56,57] and transport (slc43a1a) [58].These findings indicate that altered spawning seasons had genetic and epigenetic effects on several genes associated with cell cycle regulation.

Discussion
In the present study, we investigated the influence of altered spawning seasons on gene expression and DNA methylation in offspring.Previous studies emphasized the significant impacts of these changes on the nutritional status of both broodstock and their offspring, with no observed effects on the body weights of broodstock [4,5].Our primary aim was to reveal the genetic and epigenetic effects resulting from observed variations by utilizing liver samples from larvae at the first feeding stage.

What insights do our omics analyses provide?
Gene expression analysis revealed distinct patterns associated with altered spawning seasons.
Specifically, early and late seasons exhibited similar gene expression patterns when compared to the normal season.The differential expression of genes associated with metabolism, cellular processes, and organismal systems indicates that spawning seasons can have wide-ranging effects on biological pathways related to growth and development.
Furthermore, the analysis with overlapped DEGs revealed that all altered spawning seasons affected genes associated with metabolism, tissue development, and information processing similarly.Hence, while early and late seasons influenced a wide range of biological pathways in a similar manner, all three seasons, including the off-season, also similarly affected a certain set of genes, especially those strongly differentiated from the normal season.
DNA methylation analysis provided valuable insights into the epigenetic effects of different spawning seasons.Altered spawning seasons led to increased methylation rates, which appeared to be linked to higher body weights.Notably, these shifts primarily occurred in regions that were already highly methylated, but the rates were altered independently of the underlying methylation rates in promoter regions around transcription start sites.Our findings imply that methylation rates could serve as a measure of the current growth levels, even when comparing different treatment groups at the same developmental stage.
Several genes with multiple DMCs exhibited consistent hypo-or hyper-methylation patterns across multiple spawning seasons in promoter regions, suggesting a common epigenetic response to altered spawning seasons.These genes were associated with key biological functions, such as development and information processing, further supporting their role in influencing growth performance.Although both gene expression and DNA methylation analyses resulted in similar biological functions impacted by altered spawning seasons, they showed little overlap between DEGs and genes with DMCs.This could be attributed to regulatory differences between gene expression and DNA methylation within the same biological pathway, or it might be a limitation caused by the restriction enzyme used in RRBS.
The observed correlations between DNA methylation and gene expression were weaker than anticipated in our liver samples.Notably, CpG sites with methylation differences exceeding 30% lacked a clear link to differential gene expression.Additionally, our findings revealed inconsistent correspondence between hypo-methylation and active gene expression and between hyper-methylation and gene suppression.These results suggest that the regulation of DNA methylation is more intricate than commonly assumed as it is known that high promoter CpG methylation reduces gene expression, while low promoter CpG methylation allows for active gene expression [59].Within promoter regions, several factors, including underlying gene expression levels, transcription factor binding sites, and chromatin remodelling, likely play a role in regulating actual gene expression.

What are the significant biological differences observed in offspring?
Off-season spawning in June, facilitated through a RAS, resulted in Atlantic salmon offspring with lower larval weights compared to normal and late seasons.Offspring eggs also exhibited reduced levels of vitamin B12 and lipid classes relative to normal and late seasons.Gene expression analysis indicated significant down-regulation of genes associated with central carbon metabolism (CCM).CCM involves several biological pathways, including the citric acid cycle, and acts as a major source of energy for growth and development.Moreover, lipids are interconnected with CCM through the utilisation of acetyl-CoA as a central metabolite [60].Among the six genes with multiple DMCs from the off-season vs. normal comparison, sterol carrier protein2-like (scp-2) contained four hyper-methylated DMCs in its promoter.The scp-2 gene encodes a transfer protein that plays a key role in intracellular lipid transport [48].Overall, RAS-based off-season spawning appeared to impact lipid-mediated regulations both genetically and epigenetically.This negative impact on lipid-related mechanisms potentially accounted for impaired growth performance observed in body weight.
Offspring resulting from early season spawning in September also exhibited lower weights at the larval stage compared to those from normal and late seasons.However, the nutritional status of offspring eggs displayed higher levels of SAM/SAH, lysine, glutamine, and alanine.Gene expression analysis revealed significant down-regulation of genes associated with cell cycle regulation, while genes involved in metabolism showed both downand up-regulation.Among the five genes that exhibited differential expression and had at least one DMC in their promoters, two of them were associated with cell-cycle regulation.
Specifically, the carboxy-terminal kinesin 2 (kifc1) gene promotes mitotic spindle assembly [54], and the aurora kinase B (aurkb) gene is involved in mitotic progression [55].The kifc1 gene was hypomethylated, while the aurkb gene was hypermethylated.Consequently, the early spawning season had an impact on cell-cycle regulation both genetically and epigenetically.However, it is challenging to hypothesise that this negative impact on cellcycle regulation was linked to impaired growth performance because similar effects were observed for the late season.
Late season spawning in January showed that offspring weights were almost as close to, but slightly heavier than those from the normal season.The nutritional status of offspring eggs also exhibited similar levels of various nutrients and metabolites to the normal season, except for higher vitamin B12 and lower B-alanine.Gene expression patterns in the late season were comparable to those observed in the early season, showing significant downregulation of genes related to cellular processes, especially cell cycle regulation.Two genes, kifc1 and aurkb, showed differential expression and had at least one DMC in their promoters, and both were associated with cell-cycle regulation.Hence, the late season showed similar growth performance, nutritional status, and global DNA methylation pattern to the normal season, while gene expression pattern was similar to the early season.
Although strong similarities were observed in gene expression patterns between early and late spawning seasons, certain genes exhibited common regulation in both gene expression and DNA methylation across all three spawning seasons.Among the four DEGs identified in the off-season, choline-phosphate cytidylyltransferase B (pcyt1b) and bone morphogenetic protein 5-like (bmp5) showed concordant down-regulation across all three spawning seasons.The pcyt1b gene encodes an enzyme involved in phosphatidylcholine biosynthesis, playing essential roles in multiple metabolic pathways [25], while the bmp5 gene is involved in signalling [33].Additionally, all DMCs found in the promoter of the cytotoxic T lymphocyte-associated protein 2 beta (ctl2b) gene were hyper-methylated across all three spawning seasons.The function of ctl2b remains unknown, but it shares similarities with genes encoding cysteine proteinase, an enzyme involved in protein breakdown within cells [61].Furthermore, among genes that were both differentially expressed and had at least one DMC in their promoters, the caprin-1 gene exhibited consistent regulation, with all DMCs being hypo-methylated and its expression being suppressed.The caprin-1 gene is associated with cell cycle regulation, playing a crucial role in cellular activation or proliferation [53].This indicates that the off-season also influenced cell cycle regulation to a certain degree.
Can we link our findings with future growth performance?
Although methylation levels roughly correspond to body weights, indicating current growth levels, we can still estimate future growth performance from all kinds of observed data.
Offspring from the off-season exhibited lower lipid levels and suppressed lipid related mechanism, which could hinder their growth potential, while offspring from the early season showed suppressed cellular processes, which are also important for growth.However, the early season had the highest methylation potential (SAM/SAH ratios), suggesting superior growth potential due to methionine's role in vital growth-related biology [10,62].Therefore, even though fish from the off-season and early seasons are expected to underperform those from the normal season in later development stages, fish from the early season potentially outperforms those from the off-season at some point in later stages.
Late and normal seasons exhibited similar overall DNA methylation rates and nutritional status measures.Nonetheless, the late season displayed suppressed expression of genes involved in cellular processes, including cell cycle regulation.Furthermore, the late season exhibited the lowest SAM/SAH ratio among the four seasons.This suggests that fish from the late season may underperform those from the normal season in later development stages.
Hence, the projected future growth levels can be represented as: off-season < early < late < normal, in comparison to the current growth levels: early ≈ off-season < normal ≈ late, without considering compensation from exogenous feeding.To enhance the accuracy of these estimations, compiling tissue and developmental stage-specific DNA methylation maps would be valuable, as specific DNA methylation levels at certain loci are likely associated with the underlying growth potential.
Why are our findings important in the context of aquaculture and nutritional studies?
Our study unveiled the significant effects of altered spawning season on offspring gene expression and DNA methylation patterns in Atlantic salmon.These findings provide valuable insights into assessing current and future growth potential, which can be challenging to determine using only nutritional status and growth performance measures.Understanding the genetic and epigenetic responses to altered spawning seasons will help develop strategies to optimize growth and production in aquaculture practices.Additionally, our study can be linked to the potential impact of increasing global ocean temperatures on the spawning of aquatic vertebrates, as temperature is a key abiotic factor influencing spawning seasons.

Conclusions
The present study emphasises the substantial impact of altering spawning seasons on gene expression and DNA methylation patterns in Atlantic salmon offspring.Our observations indicate that RAS-based off-season spawning exerts an influence on lipid-mediated regulations in offspring, while sea-pen based early and late seasons affect the regulation of the offspring's cellular processes, especially cell cycle regulation.These effects potentially play a crucial role in shaping both current and future growth performance.
These findings have important implications for aquaculture and nutritional studies, as they provide valuable tools for assessing growth potential and optimizing production strategies.Future research focusing on tissue and developmental stage-specific DNA methylation maps will enhance the accuracy of growth performance estimation and provide further insights into the genetic and epigenetic responses to altered spawning conditions.

Fig. 2 .
Fig. 2. Gene expression analysis comparing altered spawning seasons to the normal season.a) PCA plots displaying clustering patterns of four spawning seasons -off-season (red), normal (green), early (blue), and late (purple) -using 20 RNA-seq samples."All mapped" plot includes 45,153 genes, while "Top 500" includes top 500 genes with high variance.b) Volcano plots showing differential expression analysis results for three pairwise comparisons: off-season vs. normal (Off-season), early vs. normal (Early), and late vs. normal (Late).Dots represent down-regulated DEGs (blue), up-regulated DEGs (red), and non-DEGs (grey).Counts of DEG and non-DEGs were provided in the legend.c) Venn diagrams illustrating overlaps of DEGs identified between two comparisons, indicating common DEGs as well as those exclusive to a single comparison.

Fig. 4 .
Fig. 4. Distribution of methylated CpG sites across genomic regions.a) Genomic region diagram illustrating three regions: flanking regions (Flank, red), promoters (P, blue), and gene bodies (GB, green).Promoter region is divided into three subregions: P250, P1K, and P5K.Gene body includes exons and introns.b) Ridge density plots showing distribution patterns of methylated CpG sites in percentage across six genomic regions.Plot includes methylation data from the present study (larvae, green) and two other studies: post-smolt (light blue) and harvest (blue) for comparison.c) Ridge density plots showing distribution patterns of methylated CpG sites for four spawning seasons: off-season (red), early (green), normal (blue), and late (purple).Plot is divided into two ranges: 0-25% and 75%-100%.

Fig. 5 .
Fig. 5. DNA methylation differences and correlations with gene expression.a) Histograms displaying distributions of methylation rate differences in six pairwise comparisons.Hypomethylated DMCs in blue, hyper-methylated DMCs in red, non-DMCs in grey.All CpG sites in plots showed significant methylation differences.b) Ridge density plots illustrating distribution patterns of methylated CpG sites across four spawning seasons: off-season (red), early (green), normal (blue), and late (purple).CpG sites were selected if at least one site identified as DMC in any of six comparisons.c) Scatter plots showing relationship between methylation rate differences and log fold changes (LFCs) of gene expression.DMCs were selected from three comparisons: off-season vs. normal (Off season), early vs. normal (Early), and late vs. normal (Late), and also selected based on two regions, P250 and P1K.Corresponding genes were linked to these DMCs chosen, and their LFCs were transformed by normal shrinkage.Data points are color-coded to indicate

1
All hypo: All three comparisons showed hypo-methylated DMCs, Hyper & Hypo: Comparisons showed both hypo-and hyper-methylated DMCs.All hyper: All three comparisons showed hyper-methylated DMCs. 2 Either hyper or hypo-methylated DMC against the normal season.The number in brackets represent the number of DMCs.*: genes selected by filtering of multiple DMCs.All 11 genes were identified by a single comparison except for slc31a1 identified in two comparisons (off-season vs. normal and late vs. normal).Most genes (9 out of 11)

Table 1 .
List of top five DEGs identified in at least two pairwise comparisons.All three comparisons showed up-regulation.2Logfoldchanges of pairwise comparison against the normal season.*:significantLFC (adjusted p-value<0.05).3Whengene symbols were unavailable, we utilised the gene symbols from either the latest NCBI version or UniProt as substitutes.

Table 2 .
List of genes with multiple DMCs in their promoter region.

Table 3 .
List of DEGs that contain DMCs in their promoter region.Either hypo-or hyper-methylated DMC against the normal season.Parentheses indicate that the site was identified as DMC, but the corresponding genes was not identified as stringent DEG.2Either up-or downregulation against the normal season.Parentheses indicate that down-or up-regulation was determined by non-