One-carbon metabolism nutrients impact the interplay between DNA methylation and gene expression in liver, enhancing protein synthesis in Atlantic salmon

ABSTRACT

Supplementation of one-carbon (1C) metabolism micronutrients, which include B-vitamins and methionine, is essential for the healthy growth and development of Atlantic salmon (Salmo salar). However, the recent shift towards non-fish meal diets in salmon aquaculture has led to the need for reassessments of recommended micronutrient levels. Despite the importance of 1C metabolism in growth performance and various cellular regulations, the molecular mechanisms affected by these dietary alterations are less understood. To investigate the molecular effect of 1C nutrients, we analysed gene expression and DNA methylation using two types of omics data: RNA sequencing (RNA-seq) and reduced-representation bisulphite sequencing (RRBS). We collected liver samples at the end of a feeding trial that lasted 220 days through the smoltification stage, where fish were fed three different levels of four key 1C nutrients: methionine, vitamin B6, B9, and B12. Our results indicate that the dosage of 1C nutrients significantly impacts genetic and epigenetic regulations in the liver of Atlantic salmon, particularly in biological pathways related to protein synthesis. The interplay between DNA methylation and gene expression in these pathways may play an important role in the mechanisms underlying growth performance affected by 1C metabolism.

KEYWORDS:

• Nutrition
• one-carbon metabolism
• B-vitamins
• methionine
• SAM
• smoltification
• omics
• DNA methylation

Introduction

One-carbon metabolism micronutrients, or simply 1C nutrients, are essential micronutrients for healthy animal development [1], including the development over the smoltification stage of Atlantic salmon (Salmo salar) [2]. While three B-vitamins, pyridoxine (vitamin B6), folate (vitamin B9), and cobalamin (vitamin B12), play a crucial role in maintaining 1C metabolism [3,4], fish need to acquire them through their diets as they cannot produce folate and vitamin B12 on their own. The main source of these vitamins for wild salmon is small pelagic fish, but the rapid development of plant-based diets in the salmon aquaculture industry has prompted reassessments of recommended levels of micronutrients in their diets [5,6]. Nonetheless, it remains challenging to determine the optimal range of 1C nutrients, such as B-vitamins and methionine, that support healthy and robust fish growth across different conditions and life stages. Moreover, there is an urgent need to explore alternative protein sources that do not rely on fishmeal or plant-based diets to enhance the sustainability of salmon aquaculture.

Methionine is a key component in 1C metabolism, providing the methyl group required for various biological functions, such as DNA methylation [7], post-translational protein modifications [8], creatine synthesis [9], carnitine synthesis [10], endogenous choline synthesis [11,12], and polyamine synthesis [13]. Before methionine can function as a methyl donor, it needs to be converted to S-adenosylmethionine (SAM) by the enzyme methionine adenosyl transferase (MAT) [14]. After donating the methyl group, SAM is converted to S-adenosylhomocysteine (SAH). As SAH is unstable and easily converted to homocysteine, which is cytotoxic to the cells [15], it needs to be either transsulfurated or remethylated to methionine. 1C metabolism is closely related to the methionine cycle and the folate cycle as they transfer and utilize one-carbon units within cells. The former involves methionine, SAM, SAH, and vitamin B12, while the latter involves vitamin B6, folate, and B12. Vitamin B6 is necessary for transsulfuration, while folate and B12 are used for the remethylation of homocysteine [16].

The liver is a key target organ for 1C metabolism and methylation reactions. Plant-based diets were found to increase lipid accumulation in the liver and intestinal tissues [17], and sub-optimal methionine level in the diets affected hepatic triglycerides (TAGs) accumulation in Atlantic salmon [18]. In addition, salmon fed low B-vitamin diets grew fattier than those fed higher B-vitamin diets [6]. Liver lipid accumulation is associated with increased metabolic stress, energy depletion, cytokine activation, and inflammation in mammals [19–21].

In the present study, we used liver samples from Atlantic salmon fed three different dietary levels of 1C nutrients, including methionine, vitamin B6, folate, and vitamin B12. Previous findings from this feeding trial showed that a medium dosage of the 1C nutrients improved growth and reduced liver size, while a high dosage had limited contribution [2]. To gain insights into the molecular mechanisms underlying these effects, we investigated two types of omics data, namely RNA-seq for gene expression and reduced-representation bisulphite sequencing (RRBS) for DNA methylation. Our goal was to gain a better understanding of the affected biological functions and identify the interactions between them by integrating these datasets.

Materials and methods

Feeding trial and sampling

The details of experimental design and sampling methods were described elsewhere [2]. In brief, triplicate groups of Atlantic salmon were fed three diets for six weeks in fresh water, through smoltification, followed by three months on-growing period in salt water. At the end of the saltwater period, liver samples were dissected for RNA and DNA extraction and HPLC analysis for SAM/SAH. The samples were immediately frozen in liquid nitrogen. Fish were anaesthetised with Tricaine (Pharmaq) before any handling. All experimental procedures and husbandry practices were conducted in compliance with the Norwegian Regulation on Animal Experimentation and European Community Directive 86/609/EEC.

Growth measurement

All fish per tank were individually measured for body weight and fork length. Following measurement, all fish were allowed to recover in aerated water. Fulton’s condition factor (K) was calculated using K = (W/L³) * 100 where W is body weight (g), and L is fork length (mm). Hepatosomatic index (HSI) was calculated as HSI (%) = (liver weight (g)/body weight (g)) * 100. More detailed methods along with the results of growth performance have been presented elsewhere [2].

SAM/SAH analysis

SAM and SAH were analysed on reverse HPLC following deproteinization in 0.4 M HClO4 as described by Wang et al. [22]. Five liver samples from each tank were pooled and homogenized for both SAM and SAH analysis.

DNA and RNA extraction

RNA was extracted using the BioRobot EZ1 and EZ1 RNA Universal Tissue kit (Qiagen) and DNase-treated with Ambion DNA-free DNA removal kit (Invitrogen, USA) according to the manufacturer’s protocols. RNA quantity and quality were assessed using NanoDrop ND-1000 Spectrophotometer (Nanodrop Technologies) and Agilent 2100 Bioanalyzer (RNA 6000 Nano LabChip kit, Agilent Technologies).

DNA isolation was performed using the DNeasy Blood & Tissue Kit (Qiagen, Cat. No. #69506) according to the manufacturer’s protocol. Liver samples were pre-treated with RNase A (provided in the Qiagen kit, 50 ng/µL, 10 min at room temperature) immediately 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, and its quantification was performed using the Qubit High Sensitivity Assay (Life Technologies #Q32854).

Atlantic salmon genome and genomic annotation

We downloaded the reference genome (ICSASG version 2) and RefSeq data (version 100) for gene annotation and coordinates from the NCBI site (https://www.ncbi.nlm.nih.gov/assembly/GCF_000233375.1). All genome regions were annotated by four genetic regions as flanking region (Flank), promoter (P), gene body (GB), and intergenic region (IGR). Flanking regions were used to investigate potential methylation effects in distal promoter, enhancer, inhibitor, and insulator sites. Gene symbols from UniProt [23] were used when gene symbols were not available at NCBI.

Library preparation of high-throughput sequencing

Liver samples were sent to the DeepSeq sequencing facility at Nord University, Bodø, Norway for RNA-seq. The library preparation was done using an NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs). The libraries were sequenced on the NextSeq500 instrument (Illumina). Detailed methods were described elsewhere [24].

Liver samples were sent to the CeMM Biomedical Sequencing Facility, Vienna, Austria for RRBS. Extracted genomic DNA were digested by MspI prior to size selection and bisulphite conversion. RRBS libraries were sequenced on Illumina HiSeq 3000/4000 instruments. Detailed methods were described elsewhere [24].

Pre-processing of high-throughput sequencing

We used the same procedures of pre-processing for RNA-seq and RRBS data as described in a previous study [24]. In short, reads were trimmed by Cutadapt [25] and Trim Galore! (Babraham Institute). Trimmed reads were aligned to the reference genome by STAR [26] with the default parameters for RNA-seq and by Bismark [27] and Bowtie 1 [28] with the default parameters for RRBS. While the mapped reads of RNA-seq were quantified by featureCounts [29], the mapped reads of RRBS were processed by Bismark for methylation calling and CpG site extraction. A total of 27 samples were divided into three diet groups (n = 9 per group) named, 1C+, 1C++, and Ctrl, defined by different 1C nutrient levels (Figure 1 & Supplementary Table S1) for both RNA-seq and RRBS. The factoextra package (https://CRAN.R-project.org/package=factoextra) was used to perform clustering analysis including PCA prior to further analysis. A variance stabilizing transformation (VST) was performed on RNA-seq counts using DESeq2 [30] prior to PCA. A non-linear clustering method, t-distributed stochastic neighbour embedding (t-SNE) [31], was additionally used for RRBS. Moreover, tank and sex effects were checked using PCA plots for RNA-seq samples (Supplementary Fig S1).

Figure 1. Overview of experimental design and nutrient levels in the feeding trial with growth performance. (a) a schematic diagram depicting different stages, pallet sizes, diets, and sampling points throughout the trial. (b) diagram showing targeted nutrients (blue) along with SAM and SAH (green) in the 1C metabolism that includes folate and methionine cycles. (c) barplots showing levels of four 1C nutrients (B6, folate, B12, and methionine) analysed for both 3 mm and 4 mm pellets for three diet groups (ctrl, 1C+, and 1C++). (d) line plot showing the growth rates in body weights (g) for three diet groups along with p-values calculated by ANOVA.

Differential gene expression analysis

Differentially expressed genes (DEGs) of three comparisons (1C+ vs Ctrl, 1C++ vs Ctrl, and 1C++ vs 1C+) were identified by DESeq2 [30] when the adjusted p-values by Benjamini-Hochberg [32] were less than 0.05. To find DEGs strongly affected by the treatment, we used log fold changes (LFCs) from pair-wise comparisons for filtering by two thresholds, |LFC| > 1.5 and |LFC| > 2. The LFC filtration was performed by using the lfcThreshold argument of the results function of DESeq2 [30]. Up- and down-regulated genes were determined based on positive and negative LFCs relative to the reference groups. The reference groups were Ctrl for both 1C+ vs Ctrl and 1C++ vs Ctrl, and 1C+ for 1C++ vs 1C + .

For the clustering analysis with density-based spatial clustering of applications with noise (DBSCAN) [33], the results of the three comparisons were combined into a pooled set of unique DEGs. The input of DBSCAN were pair-wise differences of group-wise gene expression levels calculated from the pooled set of DEGs. To calculate the average gene expression levels for each group, the trimmed mean of M-values (TMM) normalization was calculated by edgeR [34] and then scaled by the minimum and maximum values to the range between 0 and 1. The DBSCAN clustering was performed by the dbscan package (https://CRAN.R-project.org/package=dbscan) with parameters minPts = 100 and eps = 0.5.

Functional annotation with KEGG

Both over-representation analysis (ORA) and gene set enrichment analysis (GSEA) [35] on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [36] were performed by clusterProfiler [37]. The input of ORA were gene lists of DBSCAN clusters, whereas that of GSEA was shrunken LFCs calculated by the normal shrinkage method provided by DESeq2 [30], being separately calculated for three pair-wise comparisons. P-values were adjusted by Benjamini-Hochberg [32]. Enriched pathways were identified when adjusted p-values were less than 0.05 for ORA and 0.01 for GSEA. For GSEA, pathways were selected if GSEA exhibited significance in at least two of the three pair-wise comparisons (1C+ vs Ctrl, 1C++ vs Ctrl, and 1C++ vs 1C+).

Differential methylation analysis

Prior to differential methylation analysis, the number of mapped reads with less than or equal to 10 and above the 99.9^th percentile of coverage were discarded by methylKit [38]. Differentially methylated CpG sites (DMCs) of three comparisons (1C+ vs Ctrl, 1C++ vs Ctrl, and 1C++ vs 1C+) were identified by methylKit [38] when the Q-values by the sliding linear model (SLIM) method [39] were less than 0.05 and methylation differences greater than or equal to 15%. Hypo- and hyper-DMCs were determined based on positive and negative methylation differences relative to the reference groups. The reference groups were Ctrl for both 1C+ vs Ctrl and 1C++ vs Ctrl, and 1C+ for 1C++ vs 1C + .

To filter out DMCs that were strongly affected by the treatments, we used the counts of DMCs per gene in three different regions: promoter, regulatory regions within 5000 upstream from transcription start sites (TSSs), and exons. Genes were sorted by the count of DMCs in a descendant order, and the top three genes among them were selected for the further literature analysis. Genes with higher average methylation differences were selected in the case of ties.

To identify differentially methylated regions (DMRs), average methylation rates were calculated by sliding windows with size 100 and step 100 using the tileMethylCounts function of methylKit [38]. The same procedure of identifying DMCs was applied to identify DMRs after tiled regions were defined.

Linking DNA methylation with gene expression

To investigate a broad association between gene expression and DNA methylation, we divided genes into five gene groups – None, Very low, Low, Med, and High, by gene expression levels, which were defined from scaled TMM counts. Threshold ranges used to define the groups were 0 for None, (0, 0.1], (0.1, 0.5], and (0.5, 1] for Very low, Low, Med, respectively, and > 1 for High. Linking DEGs with DMCs or DMRs was done by merging NCBI gene IDs from both datasets.

Bioinformatics analysis

In-house R and Python scripts with Snakemake [40] were used to perform high-throughput sequence analysis, basic statistical analysis, and generating figures.

Results

Atlantic salmon fed with varied 1C nutrient levels showed differential growth performance

The present study examines the impact of varied 1C nutrient levels on the gene expression and DNA methylation patterns in the liver of Atlantic salmon. The feeding trial lasted 220 days through the smoltification stage (Figure 1(a)), and three diets were used: Ctrl, 1C+, and 1C++, each containing different levels of four key 1C nutrients: methionine, folate, vitamin B6, and vitamin B12 (Figure 1(b)). While the Ctrl group was fed with the recommended levels of 1C nutrients, the 1C+ and 1C++ groups were given higher levels of these nutrients (Figure 1(c) & Supplementary Table S1). The growth performance and nutrient status of the fish were previously assessed, and the results have been summarized and discussed elsewhere [2].

At the final sampling point (S4), the fish in the 1C+ group had a significantly better growth performance than those in the other two groups, with an average weight that was ~ 50 grams greater (Figure 1(d)). Other measures, including condition factor and hepatosomatic index, also indicated a healthier and better growth for the 1C+ group. In contrast, the growth performance of the 1C++ group was only slightly better than that of the Ctrl group (Supplementary Table S2).

To uncover the underlying genomic and epigenetic regulations associated with the observed differences in growth performance over smoltification, we analysed gene expression and DNA methylation of liver samples collected at the S4 sampling point (post-smolt stage). The remaining sections of the Results will focus on the analysis of these data.

Clustering analysis revealed two strong gene expression patterns that are common to the 1C-supplemented groups

We conducted an RNA-seq analysis on 27 liver samples from three groups (n = 9) to investigate the effect of 1C nutrient levels on gene expression. We identified 22,066 expressed genes by aligning read counts to the reference genome (detailed statistics for each sample are provided in Supplementary Table S3). Principal component analysis (PCA) showed a clear separation between the 1C-supplemented groups (1C+ and 1C++) and the control group (Figure 2(a)).

Figure 2. Clustering analysis of gene expression differences among three diets. (a) a PCA plot displaying the clusters of three diet groups - ctrl (blue, semi-transparent), 1C+ (yellow, semi-transparent), and 1C++ (red, semi-transparent) – using 27 RNA-seq samples including nine samples from each group. Top 1000 high variance genes were used as input data. (b) a dot plot showing the DBSCAN result in a PCA format with two identified clusters, DEG C1 (blue) and DEG C2 (red), on the pooled set of DEGs generated from three pair-wise comparisons (1C+ vs ctrl, 1C++ vs ctrl, 1C++ vs 1C+). (c) line plots showing normalized read counts as expression levels for three diet groups with the total averages for DEG C1 (blue) and DEG C2 (red) clusters.

Subsequently, we created a pooled set of differentially expressed genes (DEGs) by selecting genes that were differentially expressed in at least one of the three pair-wise comparisons (1C+ vs Ctrl, 1C++ vs Ctrl, and 1C++ vs 1C+; counts of DEGs from each comparison were provided in Supplementary Table S4). This pooled set allowed us to directly use normalized read counts as gene expression level instead of log fold changes, which is useful to identify interesting expression patterns, for example by using a robust clustering method, such as density-based spatial clustering of applications with noise (DBSCAN).

The DBSCAN clustering of the 1268 pooled DEGs showed that 74% (939 DEGs) of them were linked to two major clusters (DEG C1 and DEG C2) consisting of 533 and 409 DEGs, respectively (Figure 2(b)). The gene expression patterns of these clusters were concordant and demonstrated that all the genes in DEG C1 were down-regulated, while all the genes in DEG C2 were up-regulated in the 1C-supplemented groups (1C+ and 1C++) compared to the control group (Figure 2(c)).

The clustering analysis strongly suggests that both medium and high levels of 1C nutrients (1C+ and 1C++, respectively) affected gene expression profiles in a similar manner in the liver of Atlantic salmon, despite the significant difference in growth performance between them.

Functional analysis identified 14 biological pathways significantly affected by 1C nutrient supplementation

To examine the impact of 1C nutrient supplementation on biological pathways, we conducted functional analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database by employing two methods, namely over-representation analysis (ORA) and gene set enrichment analysis (GSEA). ORA was performed on DEGs, while GSEA was conducted on the whole gene set (see Materials and Methods for details). In total, ORA and GSEA together identified 14 KEGG pathways that can be grouped into three categories: metabolism, genetic information processing, and cellular processes (Table 1; detailed lists are provided in Supplementary Table S5 for ORA and Supplementary Tables S6-S8 for GSEA).

Table 1. Enriched KEGG pathways by ORA and GSEA.

Category	Sub-category	KEGG pathway	ORA^a	GSEA^b
Metabolism	Global and overview maps	Biosynthesis of cofactors	Down	Down (1C++)
		Biosynthesis of amino acids	Down	Down
	Energy metabolism	Oxidative phosphorylation	–	Up
	Lipid metabolism	Steroid biosynthesis	–	Down
		Linoleic acid metabolism	Up	Up (1C+)
	Amino acid metabolism	Cysteine and methionine metabolism	Down	Down (1C++)
		Arginine biosynthesis	Down	–
	Xenobiotics biodegradation and metabolism	Drug metabolism – cytochrome P450	Down	–
		Metabolism of xenobiotics by cytochrome P450	Up	–
Genetic Information Processing	Translation	Ribosome	–	Up
		Aminoacyl-tRNA biosynthesis	–	Down (1C++)
	Folding, sorting and degradation	Protein processing in endoplasmic reticulum	Up	Up (1C+)
		Protein export	Up	Up (1C+)
Cellular Processes	Cell growth and death	Cell cycle	–	Down

^aDown: enriched in C1 cluster, Up: enriched in C2 cluster, and “–”: not enriched by ORA. ^bDown and Up: enriched by both 1C+ and 1C++ groups with negative and positive NESs, respectively, but enriched only by one group when a group name (1C+ or 1C++) is followed. “–”: not enriched by GSEA.

All the six pathways supported by both ORA and GSEA showed concordant regulations regarding up- and down-regulation. Specifically, all the four pathways related to amino acid metabolism were down-regulated, whereas the remaining two pathways related to protein processing were up-regulated. Moreover, several pathways showed opposite regulations even though they were in the same subcategories. For instance, steroid biosynthesis was down-regulated, while linoleic acid metabolism was up-regulated, and both pathways are associated with lipid metabolism. In addition, two pathways related to cytochrome P450 also showed both up- and down-regulation (Table 1).

Several pathways showed both up- and down-regulation despite belonging to the same subcategories, specifically in metabolism related pathways inducing lipid metabolism and xenobiotics biodegradation by cytochrome P450. Moreover, pathways associated with biosynthesis of amino acid, particularly those linked to cysteine, methionine, and arginine showed down-regulation. In contrast, pathways related translation, such as ribosome, protein processing, and protein export, displayed up-regulation, despite all amino-acid and protein processing related pathways being potentially associated with protein synthesis.

Supplementation of 1C nutrients strongly influenced expression of genes encoding enzymes

To identify genes strongly affected by 1C nutrients, we filtered out DEGs with two threshold values of log fold changes (LFCs), namely stringent (|LFC| > 2) and semi-stringent (|LFC| > 1.5) thresholds (see Materials and Methods for details). This approach resulted in 10 and 62 genes for stringent and semi-stringent thresholds, respectively.

Among the 10 genes identified by the stringent threshold, seven belonged to the DEG C1 cluster, indicating strong down-regulation by the 1C-supplemented groups (1C+ and 1C++) compared to the control group (Table 2). Moreover, all of them were enzyme-coding genes except for two uncharacterized genes. Specifically, two paralogs of the mat2 gene, which is a rate-limiting enzyme for SAM synthesis [14], were down-regulated in a dose-dependent manner, with the 1C++ group showing more down-regulation than the 1C+ group when compared to the control group. Additionally, two paralogs of the lnx gene, which assists in the transfer of ubiquitin to target proteins [41], were also down-regulated, but the 1C+ group was more down-regulated compared to the 1C++ group.

Table 2. List of DEGs with a stringent LFC threshold (|LFC| > 2).

					LFC^5
Gene ID^a	Gene symbol^b	Gene name	C1^c	C2^d	1C+	1C++	Pattern^e
106585265	mat2	S-adenosylmethionine synthase isoform type-2-like	Y	N	−2.76	−3.05	Down
106580346	mat2	S-adenosylmethionine synthase isoform type-2-like	Y	N	−3.31	−3.85	Down
106580033	gldc	glycine dehydrogenase (decarboxylating), mitochondrial-like	Y	N	(−2.8)	−3.08	Down
106566634	sult3	cytosolic sulfotransferase 3-like	N	N	(−1.8)	−2.35	Down
106570095	kcnn2	small conductance calcium-activated potassium channel protein 2-like	Y	N	−4.40	−4.13	Down
106560268	lnx	E3 ubiquitin-protein ligase LNX-like	Y	N	−3.45	−3.09	Down
106584115	lnx1	ligand of numb-protein X 1, E3 ubiquitin protein ligase	Y	N	−3.43	(−3.1)	Down
106605369	–	uncharacterized LOC106605369	Y	N	−3.81	–	Down
106564992	gm525	uncharacterized protein C17orf67 homolog	N	N	3.99	3.84	Up
106573666	chst6	carbohydrate sulfotransferase 6-like	N	N	(3.2)	3.47	Up

^aGene ID from NCBI.^bGene symbol from UniProt. ^c,dY/N: whether the gene belongs to the corresponding cluster. ^eLog fold changes (LFCs) for 1C+ vs Ctrl and 1C++ vs Ctrl. LFCs of DEGs without a stringent LFC threshold were presented within parentheses when available. ⁶Down: down-regulation. Up: up-regulation.

As the stringent threshold resulted in no DEGs from the 1C++ vs 1C+ comparison, we also analysed the results of the semi-stringent threshold to identify genes that were potentially responsible for the limited growth performance of the 1C++ group. Among the 62 genes identified by the semi-stringent threshold, only one gene, cyp27b1 (gene ID: 106583283, gene name: 25-hydroxyvitamin D-1 alpha hydroxylase, mitochondrial-like), was from the 1C++ vs 1C+ comparison. The cyp27b1 gene, which encodes a cytochrome P450 enzyme [42], exhibited the lowest expression level for the 1C++ group in an order of 1C++ < Ctrl < 1C+.

To summarize, 1C nutrient supplementation strongly influenced the regulation of the genes that encode various enzymes, including SAM synthesis ubiquitination, and cytochrome P450. Moreover, altered gene expression of the mat2 gene might have affected the regulation of SAM synthesis, which potentially influenced DNA methylation profiles even though the direct association between SAM synthesis and DNA methylation profiles is poorly understood.

The SAM/SAH ratio indicated a low methylation capacity in fish given a high dosage of 1C nutrients

The SAM/SAH ratio is commonly used as marker to assess methylation capacity [43], and both SAM and SAH play a crucial role in the methionine cycle (Figure 1(b)). We performed reversed-phase chromatography to measure SAM and SAH levels in three feeding groups in the liver of Atlantic salmon and found that the SAM level was significantly lower in the 1C++ group compared to the other two groups. Although SAH level was higher in the 1C++ group, the difference was not statistically significant. The SAM/SAH ratio for the 1C++ group was significantly lower than that for the other groups, indicating lower methylation capacity for the 1C++ group (Figure 3(a)).

Figure 3. Global and regional DNA methylation landscapes along with SAM/SAH ratio. (a) box plots showing SAM and SAH levels (nmol/g) with SAM/SAH ratio for three diet groups. P-values are calculated by ANOVA, and red stars on the top indicate statistical significance from the base mean by t-test. Dotted horizontal lines indicate base means. (b) a t-SNE plot displaying the clusters of three diet groups - ctrl (blue, semi-transparent), 1C+ (yellow, semi-transparent), and 1C++ (red, semi-transparent) – using 27 RRBS samples. (c) a ridge density plot showing average methylation rates of 157,201 CpG sites for three diet groups. (d) a schematic diagram showing the definition of three different genomic regions – gene body (GB, green), promoter (P, blue), and flanking regions (flank, red) along with intergenic regions (IGRs). (e) ridge density plots showing average methylation rates for three diet groups in four different genetic regions. (f) a running average line plot showing average DNA methylation rates of all RRBS samples within 5000 bp up- and down-stream around TSSs. The running average lines were calculated separately for five different genetic regions – RS (red), P (dark yellow and green), and GB (blue and pink).

These findings suggest that a high dosage of 1C nutrients (1C++) may impair methylation capacity, which can have implications for the lower growth performance of the 1C++ group compared to the 1C+ group.

Supplementation of 1C nutrients broadly contributed to the shifting of DNA methylation rates that are positively correlated with growth performance

We conducted an RRBS analysis on liver samples from three groups (n = 9) to examine the effect of 1C nutrient levels on DNA methylation profiles. After aligning the sequence reads to the salmon genome (see Supplementary Table S9 for alignment statistics), we identified over 150,000 CpG sites that were methylated.

Despite using various clustering approaches, including t-distributed stochastic neighbour embedding (t-SNE), we found no clear clustering of the CpG sites among the three feeding groups (Figure 3(b) for t-SNE and Supplementary Fig S2 for other clustering results). This implies that the individual variability in methylation profiles was higher than the group variability.

However, we observed noticeable peak shifts in methylation rates (Figure 3(c)), with the distributions of methylation rates being significantly shifted in an order of Ctrl < 1C++ < 1C+ (Kolmogorov – Smirnov test; Supplementary Table S10). This shift was consistent with the observed growth performance (e.g., mean body weight shown in Figure 1(d)). Additionally, we found that the shift was more pronounced for CpG sites with high methylation rates (Figure 3(c), between 75% and 100%), indicating that these highly methylated sites were even more methylated in the 1C-supplemented groups (1C+ and 1C++) compared to the control group (Figure 3(c)).

Shifting patterns of DNA methylation rates were similar in different genomic regions

To investigate region-specific methylation patterns, we annotated CpG sites into four genomic regions (Figure 3(d)): gene body (GB), promoter (P), flanking regions (Flank), and intergenic region (IGR). Similar to the shift that we observed at the whole genome level, the distributions of methylation rates were significantly shifted in an order of Ctrl < 1C++ < 1C+ in all four genetic regions, except for the comparison between 1C+ and 1C++ in the promoter region that showed no significant difference (Figure 3(e) and Supplementary Table S10).

The underlying methylation rates appeared to be different in the promoter regions compared to the other regions, as more than half of the CpG sites were associated with low methylation rates (Figure 3(e), between 0% and 25%). This can be explained by low methylation rates around transcription start sites (TSSs). We observed that methylation rates decreased sharply to around 25%, with the lowest point at the transcription start sites (TSSs), from an average methylation rate of approximately 80% in the surrounding area (Figure 3(f)).

Our findings suggest that 1C nutrient supplementation modulated methylation profiles similarly in different regions but with an exception in the regions around TSSs, which showed low methylation rates and less shifting.

Differential methylation analysis revealed more hyper-methylated CpGs in the 1C-supplemented groups when compared to the control group

We performed a differential methylation analysis to identify CpG sites that were differentially methylated with statistical significance between treatment groups. We used three pair-wise comparisons (1C+ vs Ctrl, 1C++ vs Ctrl, and 1C++ vs 1C+), which lead to over 6000 differentially methylated CpG sites (DMCs) for each comparison (Supplementary Table S11).

The result showed that approximately 70% of DMCs identified in two comparisons (1C+ vs Ctrl and 1C++ vs Ctrl) were hyper-methylated (Figure 4(a)). Additionally, the 1C+ group had slightly more hyper-methylated DMCs than the 1C++ group (shown as hypo-methylated DMCs in Figure 4(a) as the 1C+ group was used as a reference/base group). This finding is consistent with our observation that methylation rates were shifted in the order of Ctrl < 1C++ < 1C+ with the mean methylation rates of 83%, 84.2%, and 84.5%, accordingly.

Figure 4. Distributions and counts of DMCs identified by three pair-wise comparisons. (a) violin plots showing distributions and counts of non-DMCs, hypo DMCs, and hyper DMRs in four different genomic regions for three comparisons. The x-axis represents the percentage differences of methylation rates. (b) stacked barplots showing the number of DMCs per gene, 1 DMC (red), 2 DMCs (green), and greater than or equal to 3 DMCs (blue), in three different genomic regions for three comparisons. The x-axis represents the proportions of DMC counts in percentage.

Filtering multiple DMCs identified key genes affected by 1C supplementation

Our analysis revealed that most genes containing DMCs had only a single DMC, with more than 55%, 68%, and 66% of genes supported by single DMCs in gene bodies, promoters, and regulatory sequences, respectively (Figure 4(b)). While a single DMC may be sufficient to regulate gene expression [44], it could also be a false positive or a single nucleotide variant. Therefore, we filtered out genes that contained multiple DMCs, based on several criteria, and selected the top three genes from six different categories (see Supplementary Tables S12-S14, and Materials and Methods for details).

Our filtration resulted in 22 genes containing multiple DMCs in their genetic regions: promoters (P), 1K-5K upstream from the transcription start site (RS (5K)), and exons (GB) (Table 3). All genes except one (tbx15) contained either hypo-methylated or hyper-methylated DMCs but not both. More than a quarter of the genes (6 out of 22) were identified by both the 1C+ vs Ctrl and 1C++ vs Ctrl comparisons, indicating that 1C supplementation altered the DNA methylation of these six genes (rn182, smu1-like, brf1, rnf186, c3, and one uncharacterized gene) similarly regardless of the dosage levels (Table 3).

Table 3. List of top three DMCs identified by three different comparisons.

Region^a	Gene symbol^b	Gene name	1C+ vs Ctrl^c	1C++ vs Ctrl^c	1C++ vs 1C+^c	Function
P	rn182	RING finger protein 182	6 (hyper)	4 (hyper)		Ubiquitination
	-	uncharacterized LOC106572013	4 (hyper)		4 (hypo)	-
	styxl1	serine/threonine/tyrosine-interacting protein A-like	3 (hyper)			Translation, Signalling
	loxl3	lysyl oxidase homolog 3-like		3 (hyper)
	yars	tyrosyl-tRNA synthetase		4 (hypo)		Translation
	top2b	DNA topoisomerase 2-beta-like			2 (hyper)	Immune
	nca11	Neural cell adhesion molecule 1-A			3 (hypo)	Cell adhesion
Flank (5K)	vngl	vang-like protein 1	5 (hyper)			Signalling
	-	uncharacterized LOC106577634	6 (hyper)	4 (hyper)		-
	smu1-like	WD40 repeat-containing protein SMU1	1 (hyper)	3 (hyper)	2 (hyper)	Transcription
	hgf	hepatocyte growth factor-like		1 (hyper)	4 (hyper)	Signalling
	csrnp3	cysteine/serine-rich nuclear protein 3-like	5 (hypo)		2 (hyper)	Transcription
	brf1	transcription factor IIIB 90 kDa subunit-like	2 (hypo)	4 (hypo)		Transcription
	btg4	B-cell translocation gene 4	2 (hypo)		8 (hyper)	Cell cycle
	tbx15	T-box 15			3 (both)	Adipose tissue regulation
GB (exon)	rnf186	RING finger protein 186-like	6 (hyper)	4 (hyper)		Ubiquitination
	gcnt3	beta-1,3-galactosyl-O-glycosyl-glycoprotein beta-1,6-N-acetylglucosaminyltransferase 3-like		3 (hyper)	4 (hyper)	Cell adhesion
	st8sia4	ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 4	5 (hypo)		6 (hyper)	Cell adhesion
	xirp2	xin actin-binding repeat-containing protein 2-like	5 (hypo)			Cell adhesion
	cd276	CD276 antigen-like		5 (hypo)	4 (hypo)	Immune
	c3	complement C3-like	4 (hypo)	4 (hypo)		Immune
	lbh	protein LBH-like	4 (hyper)	2 (hypo)	4 (hypo)	Signalling

^aP: promoter regions, Flank (5K): Flaking regions between 1K and 5K upstream from TSS, GB (exon): exons. ^bGene symbols from UniProt. ^cNumber of DMCs identified by three pair-wise comparisons. hypo, hyper, both: all DMCs are hypo- or hyper-methylated or both hypo and hyper-methylated, respectively.

Literature analysis on these 22 genes revealed that they were associated with eight biological functions: ubiquitination (rn182 and rnf186) [45], translation (styxl1 and yars) [46,47], transcription (smu1-like, csrnp3, and brf1) [36,48,49], signalling (styxl1, vngl, hgf, and lbh) [36,46], immune system (top2b, cd276, and c3) [50–52], cell cycle (btg4) [53], cell adhesion (nca11, gcnt3, st8sia4, and xirp2) [54–56], and adipose tissue regulation (tbx15) [57].

It is important to note that none of these 22 genes appeared to be differentially expressed at the S4 sampling point (post-smolt stage). These results suggest that 1C nutrients affected DNA methylation profiles of genes involved in various cellular processes without directly changing their gene expression levels in the liver.

Integration of omics data revealed distinct DNA methylation landscapes in promoters and gene bodies depending on gene expression levels

To investigate the overall relationship between DNA methylation and gene expression, we integrated RNA-seq and RRBS data using all available genes and CpG sites. Most genes showed weak expression levels in the liver, as indicated by the right-skewed distribution of gene expression (Figure 5(a)). By considering this skewness, we grouped genes into five sets based on their normalized read counts (Figure 5(b)) by using four threshold values (shown as red dotted lines in Figure 5(a)). We also calculated average methylation rates from all 27 RRBS samples and defined six genetic regions: two each from promoters (P), flanking regions (Flank), and gene bodies (GB).

Figure 5. DNA methylation landscapes in five gene expression groups and six different genomic regions. (a) a density plot showing the distribution of gene expression calculated from normalized read counts from RNA-seq samples. Red vertical lines represent threshold values to define five different gene groups. (b) a barplot showing the count of genes in the five gene expression groups - none, very low, low, med, and high – defined by different gene expression levels. The ‘none’ group contains genes without any expression. (c) ridge density plots displaying distributions of methylation rates for five different gene groups in six genetic regions.

Our results show that promoters and gene bodies had distinct DNA methylation landscapes depending on the expression levels of associated genes (Figure 5(c)). Promoter regions of higher expression genes showed lower methylation rates (0–25%), especially when located near the transcription start site (TSS). Conversely, promoter regions of genes without expression exhibited both high and low methylation rates. This trend was more pronounced when the regions were located closer to the TSS (as shown by comparing P (250) and P (1K) in Figure 5(c)). Furthermore, we observed that highly methylated regions (75–100%) were even more methylated (as peaks are shifted towards 100% in Figure 5(c)) in the bodies of genes with stronger expression levels, with a more pronounced trend observed in introns than in exons.

Overall, our findings suggest that DNA methylation landscapes vary depending on the location of the genome and the expression levels of associated genes, particularly in promoters and gene bodies. These results provide insight into the complex interplay between DNA methylation and the regulation of gene expression regulation.

Both hypo- and hyper-methylated regions showed complex relationships with gene expression in promoters and exons

To specifically identify genes with significant differences in gene expression and DNA methylation induced by 1C nutrients, we calculated differentially methylated regions (DMRs) and selected those located within 1000 up and downstream from TSSs. This approach identified nine genes that were differentially expressed with differentially methylated regions around TSS (Table 4, Supplementary Table S15).

Table 4. List of genes that have DMRs around their TSSs.

DMR^a	Gene symbol^b	Gene name	Treat^c	Region^d	Rate range^e	DEG^f	Expr^g	Function
Hypo	hhipl1	HHIP-like protein 1	1C++	P	Low	Down	High	Metabolism
	slc51a	organic solute transporter subunit alpha-like	1C+	P	Low	Down	High	Transport
	fyb1	FYN-binding protein 1	1C+	P	Low	Up	Low	Signalling
Hyper	pla2g10	group 10 secretory phospholipase A2-like	1C+	P	High	Up	Low	Metabolism
	prelid1	PRELI domain-containing protein 1, mitochondrial-like	1C+, 1C++	P	High	Up	High	Apoptosis
	mpst	3-mercaptopyruvate sulfurtransferase-like	1C+	Exon	Low	Up	High	Metabolism
	angptl3	angiopoietin-like 3	1C+	Exon	Med	Down	Med	Cell adhesion
	slc3a1	solute carrier family 3 (amino acid transporter heavy chain), member 1	1C++	Exon	High	Down	Med	Transport
	cyp2m1	cytochrome P450 2M1-like	1C++	Exon	High	Down	Med	Metabolism

^aHypo: Hypo-methylated DMR, Hyper: hyper-methylated DMR. ^bGene symbols from UniProt. ^cTreatment group compared to the control group used for the DMR identification. ^dP: promoter region, and Exon: exon region within 1000 bp from TSS. ^eThree categories of methylation shift indicating that the shift occurred within a range of 0%-50% for Low, 50%-100% for High, and 25%-75% for Med. ^fDown: downregulated DEGs, and Up: up-regulated DEGs against the control group. ^gGroups formed by different gene expression levels: High, Med, and Low.

These genes were both hypo- and hyper-methylated in their promoters, but only hyper-methylated in the exons. In promoters, the methylation shifts occurred in a lower range of methylation rates (<50% or shown as ‘Low’ in the Rate range column in Table 4) for the hypo-methylated DMRs, while they occurred in a higher range (>50%) for the hyper-methylated DMRs. However, in exons, the shifts of methylation occurred in a wide range of methylation rates. In addition, gene expression levels varied among the nine genes (the Expr column in Table 4). Notably, these nine genes displayed no common patterns in terms of positive or negative correlation between gene expression and DNA methylation in the liver.

In addition, literature analysis on the nine genes revealed that they were associated with five biological functions: metabolism (hhipl1, pla2g10, mpst, and cyp2m1) [36,58,59], cell adhesion (angptl3) [60], signalling (fyb1) [61], transport (slc51a and slc3a1) [62,63], and apoptosis (prelid1) [64].

These DMRs were potentially involved in several different categories of epigenetic regulations, suggesting that they may have a complex role in gene expression rather than simply functioning as canonical on-off switches of DNA methylation in animals.

Discussion

The growth performance of Atlantic salmon was highest with the medium dosage of 1C nutrients (1C+), followed by the high dosage (1C++) and the control diet [2]. In a similar Atlantic salmon study utilizing a feeding trial with three different levels of 24 micronutrient components, including B-vitamins (B6, B9, and B12), the growth performance was different as high and medium > low dosage [24]. The study also exhibited that the micronutrients affected transcriptional and epigenetic patterns in a dose-dependent manner. The duration of this study was approximately 54 weeks compared to around 31 weeks of our study. Interestingly, in this 24-micronutrient study, the growth performance was higher in the medium dosage group compared to the high and low dosage groups at around the 10th week. Hence, if our study had also continued further, there could be a possibility that the growth performance of the 1C++ group could outperform that of the control group.

Both 1C supplementation induced significant changes in gene expression across a wide range of biological pathways in a similar manner in the liver. Specifically, the supplements affected pathways related to metabolism, genetic information processing, and cellular processes (Table 1). However, up- and down-regulation of the affected genes varied within related pathways. For instance, while both supplements led to down-regulation of genes involved in steroid biosynthesis, they induced up-regulation of genes related to linoleic acid metabolism, even though both pathways are part of the lipid metabolism. Similarly, inconsistencies were found in pathways related to cytochrome P450 and translation. Moreover, the 1C supplements had an interesting effect on the amino-acid biosynthesis pathways and protein synthesis. Specifically, the supplements appeared to reduce the biosynthesis of certain amino acids while increasing protein synthesis. Furthermore, the study summarizing the growth and nutritional status of the same experiment as the current study revealed average protein productive values (PPV, %) ± SEM of 46.4 ± 1.2, 51.1 ± 0.1, and 52.2 ± 2 for the Ctrl, 1C+, and 1C++ groups, respectively [2]. These values, used to estimate protein production, suggest that the 1C supplement groups enhanced protein synthesis even though the difference was not statistically significant. This implies that the liver tissues may enhance protein synthesis without needing to activate specific amino-acid biosynthetic pathways when 1C supplements were provided.

Although most genes significantly affected by the 1C supplements showed similar expression patterns for both 1C-supplemented groups, certain genes, such as mat2, lnx, and cyp27b1, exhibited divergent patterns (see Table 2 for mat2 and lnx homologues). The 1C supplements reduced the expression of two mat2 paralogs, which were more down-regulated in the 1C++ group than the 1C+ group. Since mat2 encodes a rate-limiting enzyme for SAM synthesis [14], the down-regulation of mat2 could be associated with the low methylation capacity of the 1C++ group, as indicated by the SAM/SAH ratio. The 1C supplements also reduced the expression of two lnx paralogs, with the 1C+ group being more down-regulated than the 1C++ group. The lnx gene encodes an E3 ubiquitin ligase, which plays a key role in transferring ubiquitin to various protein substrates [41]. Ubiquitination is a tightly regulated post-translational modification that marks proteins for a wide range of functions, including protein degradation and subcellular localization. Furthermore, n-terminal methionine-linked ubiquitin, or linear ubiquitin, is a significant post-translational modification involved in innate immune signalling [65]. Lastly, the cyp27b1 gene, which encodes a cytochrome P450 enzyme, exhibited the highest expression level in the 1C++ group, followed by the control and 1C+ groups. This enzyme plays an important role in calcium homoeostasis by catalysing the active form of vitamin D [42]. Moreover, many other enzymes belonging to the P450 superfamily are related to cholesterol and steroid regulations [66]. In summary, these three genes and their associated regulations, such as the methylation capacity, ubiquitination, and calcium homoeostasis, are potential contributors to the observed growth performance in our study even though the expression patterns of these genes are not completely in line with the growth performance pattern.

While most DNA methylation markers are established during early developmental stages and become stable afterwards, some markers are dynamic throughout the entire life stage [7,67]. Our DNA methylation analysis revealed that both hypo- and hyper-methylation markers appeared to be well-conserved, as we observed very few methylation markers shifting from low methylation (0–25%) to high methylation (75%-100%) or vice versa (Figure 4(a)). However, we also noticed that already methylated sites became even more methylated with 1C supplementation in the same order as the observed growth performance (Figure 3(c,e) and Supplementary Table S10). This suggests that these shifted markers can be functional as DNA methylation markers for assessing growth performance.

Gene annotation of DMCs and filtering with DMC counts led us to identify 22 genes with multiple DMCs (Table 3). Functional annotation of these genes revealed that they were associated with seven different biological categories: ubiquitination, translation, transcription, signalling, immune system, cell adhesion, and adipose tissue regulation. Among these categories, we focused on ubiquitination and translation, as their gene expression levels were also affected by the 1C supplements. For genes associated with ubiquitination, we observed hyper-methylation of two RING finger genes (rn182 and rnf186) and significant down-regulation of gene expression in two other genes (lnx and lnx1). RING finger genes and lnx homologues are all E3 ubiquitin ligases linked to ubiquitination [41,45]. For genes associated with translation, we identified two enzyme-encoding genes (styxl1 and yars). The styxl1 gene plays crucial roles in several cellular functions, including protein synthesis, signal transduction, and cell growth [46], while the yars gene encodes tyrosyl-tRNA synthetase that catalyses the aminoacylation of tRNA [47]. Notably, our KEGG analysis revealed that the expression levels of genes involved in the aminoacylation of tRNA were significantly down-regulated in the 1C++ group (Table 1), and the promoter of the yars gene was hypo-methylated in the same group (Table 3). These findings suggest that the 1C supplements influenced both gene expression and DNA methylation status of the genes involved in protein synthesis and ubiquitination.

Instead of limiting our analysis to genes showing significant differences, such as DEGs and genes with DMCs or DMRs, we took a more comprehensive approach by studying DNA methylation across various genomic regions and gene expression levels, by utilizing the whole datasets (Figure 5). Our findings highlighted significant differences in methylation patterns, particularly in promoters, and some variations in gene bodies. Specifically, we observed that genes with higher expression levels exhibit lower methylation rates near their transcription start sites, while their gene bodies tend to have increased methylation in regions that were already highly methylated. These findings have raised interesting questions about the relationship between DNA methylation and gene expression. While it is commonly believed that methylation in promoters is linked to gene repression and that hypo-methylation in promoters and hyper-methylation in gene bodies are associated with gene expression [7], these correlations may not always be valid across different regions and gene expression levels. It is important to note that this hypothesis may have a strong impact on how the results of nutritional experiments should be interpreted as differences among treatment groups of such experiments can be subtle. Therefore, it is crucial to provide additional information when combining gene expression and DNA methylation data instead of solely relying on correlation coefficients for filtering.

Prior to combining two omics data sets, we implemented two pre-processing steps: identification of DMRs and filtering of DMRs to only include those around TSSs. This reduced noise and restricted the number of candidate genes whose gene expression and DNA methylation were significantly influenced by 1C supplements. By merging the two omics data sets, we identified nine genes associated with various biological functions, such as metabolism, signalling, cell adhesion, and transport (Table 4). We did not observe any consistent correlation patterns between gene expression and DNA methylation for these genes, suggesting their potential involvement in complex transcriptional regulation through enhancers, silencers, and chromatin remodelling. Furthermore, some of these genes were implicated in multiple biological pathways. For example, pla2g10 and mpst, two enzymes linked to metabolism, were associated with a wide range of KEGG pathways, including cysteine and methionine metabolism, and linoleic acid metabolism [36]. The slc51a gene is responsible for bile acid transport but also plays a role in steroid transport [62], while the cyp2m1 gene encodes a cytochrome P450 enzyme with fatty acid hydroxylation activity [58].

Our study shares noticeable similarities with a previous salmon study that investigated the effects of varying levels of 24 micronutrients on growth performance [24]. Despite differences in feeding trials, both studies found a positive association between the down-regulation of steroid synthesis in the liver and improved growth. Additionally, we observed dynamic methylation alterations in CpG sites associated with genes involved in cell signalling and cell adhesion (Table 3), which is consistent with the previous study’s findings. Furthermore, our results align with another study that examined metabolites and gene expression in skeletal muscle tissues from the same feeding trial of 1C nutrients [68]. This study demonstrated that the down-regulation of amino acid biosynthesis was associated with reduced levels of free amino acids in the muscle tissues of the 1C-supplemented group.

The effects of 1C nutrients go beyond their role in fatty acid metabolism and cell adhesion and signalling. Specifically, we found that the groups supplemented with 1C nutrients potentially had the ability to bypass the need for synthesizing certain essential amino acids, such as methionine, and allocate more resources towards protein synthesis (Table 1). As unnecessary and damaged proteins are tagged by ubiquitin for degradation [69], the control group potentially required recycling of proteins more than the 1C-supplemented groups. Additionally, a high dosage (1C++) had a negative impact on methylation capacity (from SAM/SAH in Figure 3(a)). Although the underlying mechanism is unclear, the difference of DNA methylation patterns between two 1C-supplemented groups may have been associated with better growth performance of the medium dosage group. Furthermore, only a few genes showed significant differences in expression and DNA methylation at the same time, many affected genes were associated with the same underlying biological pathways and regulations.

Conclusions

Our findings provide valuable insights into the complex regulation of gene expression and DNA methylation that are associated with growth performance, particularly with regard to protein synthesis. These insights could have significant implications for developing strategies to enhance growth and improve health outcomes from both genetic and epigenetic perspectives in salmon aquaculture, as well as in a wide range of nutrient studies.

Author contributions

K.H.S conceived and designed the research. M.E., V.V. and T.H.T. conducted and designed the feeding trial. J.M.O.F. and C.B. prepared RNA-seq and RRBS data. T.S., M.E., A.A., and K.H.S. analysed and interpreted data. T.S. drafted the manuscript, and J.M.O.F, A.A., and K.H.S substantially revised it. The authors read and approved the final manuscript.

Acknowledgments

We thank Amelie Nemc and Bekir Ergüner for advice on RRBS analysis, and Hui-Shan Tung and Eva Mykkeltvedt for DNA and RNA extraction at the Institute of Marine Research. This work was supported by The Research Council of Norway under project no: 267787 (NutrEpi) and by the Institute of Marine Research under the Nutritional programming project.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Raw read data for both RNA-seq samples (PRJNA680426) and RRBS samples (PRJNA642024) have been uploaded to the Sequence Read Archive (SRA) on the NCBI website. To ensure easy access to these resources, we have prepared seven data files in a tabular format, which have been uploaded to a public data repository (https://figshare.com/s/9ab1ed68cfa1debe29d4). Additionally, basic descriptions of the data and columns can be found in Supplementary Tables S16-S22.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2024.2318517

Additional information

Funding

The work was supported by the Havforskningsinstituttet [Nutritional programming]; Norges Forskningsråd [267787].