Skip to main content
Download PDF

Broodstock nutritional programming differentially affects the hepatic transcriptome and genome-wide DNA methylome of farmed gilthead sea bream (Sparus aurata) depending on genetic background

BMC Genomics volume 24, Article number: 670 (2023) Cite this article

Abstract

Background

Broodstock nutritional programming improves the offspring utilization of plant-based diets in gilthead sea bream through changes in hepatic metabolism. Attention was initially focused on fatty acid desaturases, but it can involve a wide range of processes that remain largely unexplored. How all this can be driven by a different genetic background is hardly underlined, and the present study aimed to assess how broodstock nutrition affects differentially the transcriptome and genome-wide DNA methylome of reference and genetically selected fish within the PROGENSA® selection program.

Results

After the stimulus phase with a low fish oil diet, two offspring subsets of each genetic background received a control or a FUTURE-based diet. This highlighted a different hepatic transcriptome (RNA-seq) and genome-wide DNA methylation (MBD-seq) pattern depending on the genetic background. The number of differentially expressed transcripts following the challenge phase varied from 323 in reference fish to 2,009 in genetically selected fish. The number of discriminant transcripts, and associated enriched functions, were also markedly higher in selected fish. Moreover, correlation analysis depicted a hyper-methylated and down-regulated gene expression state in selected fish with the FUTURE diet, whereas the opposite pattern appeared in reference fish. After filtering for highly represented functions in selected fish, 115 epigenetic markers were retrieved in this group. Among them, lipid metabolism genes (23) were the most reactive following ordering by fold-change in expression, rendering a final list of 10 top markers with a key role on hepatic lipogenesis and fatty acid metabolism (cd36, pitpna, cidea, fasn, g6pd, lipt1, scd1a, acsbg2, acsl14, acsbg2).

Conclusions

Gene expression profiles and methylation signatures were dependent on genetic background in our experimental model. Such assumption affected the magnitude, but also the type and direction of change. Thus, the resulting epigenetic clock of reference fish might depict an older phenotype with a lower methylation for the epigenetically responsive genes with a negative methylation-expression pattern. Therefore, epigenetic markers will be specific of each genetic lineage, serving the broodstock programming in our selected fish to prevent and mitigate later in life the risk of hepatic steatosis through changes in hepatic lipogenesis and fatty acid metabolism.

Peer Review reports

Background

Many studies in humans and animal models have demonstrated that sub-optimal nutrition during pregnancy and neonatal stages induced metabolic changes that can manifest at the tissue, cellular and molecular levels, leading marked physiological consequences for the offspring [1, 2]. Indeed, nutritional stresses act on genes or gene pathways common to most insults (gatekeeper genes), and such knowledge has contributed to prevent diseases in humans [3] or improve performance in sheep [4] or beef cattle [5]. In fish, early nutritional programming also results in developmental adaptations [6,7,8,9,10], including the improved acceptance and utilization of plant-based diets in a typically marine fish such as gilthead sea bream [11,12,13]. Such nutritional intervention operates among other paths through changes in the offspring hepatic lipid metabolism, being attention initially focused on fatty acid desaturase 2 (Fads2) that shared a clear functional diversification across fish species [14, 15]. Certainly, Fads2 (Δ6-desaturase) catalyses the first and rate limiting step in the biosynthesis of n-3 long-chain polyunsaturated fatty acids (LC-PUFA) to convert α-linolenic acid (ALA, 18:3n-3) into eicosapentaenoic acid (EPA, 20:5n-3). However, the loss of Fads1 (Δ5-desaturase) in marine fish blocks the availability of these animals to elongate and desaturate EPA until docosahexaenoic acid (DHA, 22:6n-3). Alternatively, Fads2 of some marine herbivorous fish is able to produce DHA from EPA through the Δ4 desaturation pathway [16,17,18,19]. Such discovery highlighted that the fish biosynthetic capacity of n-3 LC-PUFA not only depends on the dichotomy between freshwater and marine fish species, but also on the trophic level [20].

Altogether, the above findings reinforce the role of Fads2 as a rate-limiting step in the biosynthesis of n-3 LC-PUFA in marine fish [21, 22], and selective breeding for enhanced broodstock fads2 expression improved the offspring utilization of plant-based diets (limited supply of n3-LC-PUFA) in gilthead sea bream (Sparus aurata) [23]. Otherwise, de novo fatty acid biosynthesis (Δ9-pathway) offers the possibility to mitigate the signs of deficiencies in essential fatty acids through the increased production of monounsaturated fatty acids (oleic acid, 18:1n-9) instead of EPA and DHA. Such adaptive feature will contribute to preserve the fatty acid unsaturation index of membrane phospholipids, and juveniles of gilthead sea bream fed semi-synthetic diets formulated to be deficient in n-3 LC-PUFA shared a marked up-regulated expression of the novo hepatic lipogenic genes, elongase 6 (elovl6) and stearoyl-coenzyme A desaturase 1a (scd1a) [24]. However, the risk of hepatic steatosis cannot be underestimated with exaggerated or poorly regulated hepatic lipogenesis [25], and broodstock nutritional programming with an enriched-ALA diet served to maintain regulated the enhanced expression of scd1a in the gilthead sea bream offspring. This thing resulted in a negative correlation between the hepatic scd1a expression and the DNA methylation level of several CpG sites of a CG island of the proximal promoter region that contains a PUFA responsive element [24]. All this supports the highly epigenetic regulated expression of scd1a in gilthead sea bream, though further research is needed when considering the extent to which such response might involve other processes, directly or indirectly related to lipid metabolism, and more importantly how all this can be driven by a different genetic background.

The aim of this study is to contribute to solve the gap of knowledge for the interplay between epigenetics and genetics in gilthead sea bream, combining massive gene expression analysis (RNA-seq) with genome-wide DNA methylation approaches, for which the gold-standard technique is the whole genome bisulphite sequencing (WGBS). This precise method provides a single-base CpG resolution [26,27,28], but sequencing to a sufficient coverage can result economically unaffordable. Reduced representation bisulphite sequencing is a cost-effective alternative [29, 30] that has been used in aquaculture for targeting the CG-rich genomic regions of a wide range of species, including European sea bass [31, 32], Atlantic salmon [33], rainbow trout [34], and Nile tilapia [35]. However, bisulphite conversion-based assays fail to differentiate between 5-methylcytosine and other epigenetic modifications [36,37,38]. Instead, despite not allowing for a single-base CpG resolution, the methyl-binding domain sequencing (MBD-seq) offers an economical alternative for a large coverage of the CpG methylome [39,40,41] that approximates the sensitivity/specificity of WGBS [42]. Such approach was successfully employed in human and animal studies [43,44,45,46], including fish [47, 48] and it was applied herein to underscore how broodstock nutrition with low fish oil feed formulations was affected by selective breeding within the PROGENSA® gilthead sea bream program, which selected for fast growth [49] and a low incidence of skeletal deformities [50], but also for changes in behaviour and swimming performance [51, 52], and intestinal microbiota plasticity to cope with changes in diet composition [53,54,55].

Results

Fish performance

Data of fish performance are shown in Table 1. Initial body weight was similar in all animals, regardless of diet and genotype. Reference (REF) fish presented a significantly lower final body weight (BW) compared to genetically selected (GS) fish (P < 0.001). Similarly, diet affected gilthead sea bream weight gain, presenting fish fed FUTURE diet lower final BW than fish fed control (CTRL) diet (P < 0.001) within the same genotype. At the end of the experimental period, as shown by the daily growth index (DGI), GS fish fed the CTRL diet grew significantly (P < 0.05) better than the animals used for the rest of the treatments, whereas GS fish fed the FUTURE diet performed similarly to those REF fish fed the CTRL diet, and REF fish fed the FUTURE diet presented the lowest (P < 0.05) values. Statistical analysis revealed also a significant effect of genotype and diet on specific growth rate (SGR; P < 0.001). Additionally, either fish genotype or diet formula affected diet utilization, presenting GS gilthead sea bream better feed conversion ratio (FCR) than the REF animals (P < 0.001) as well as fish fed the CTRL diet showed lower FCR than fish fed FUTURE diet (P < 0.05).

Table 1 Growth and feed utilization parameters in genetically selected fish (GS) and reference fish (REF) fed either FUTURE or CTRL diets during the challenge phase
Full size table

Patterns of offspring gene expression

Illumina sequencing of offspring mRNA liver samples from crosses of REF or GS fish challenged with a CTRL or FUTURE diet generated ~ 2,008 million paired-end (PE) reads (2 × 150), with an average of ~ 83 million reads per sample (Additional file 1: Supplementary Table 1). After trimming and quality filtering, around 9% of all liver reads were discarded, and the remaining reads ranged between 49 million (~ 7.35 Gb) and 106 million (~ 15.9 Gb) in all samples. Up to 92% of these pre-processed reads were mapped against the IATS-CSIC gilthead sea bream reference genome, which retrieved 46,545 expressed coding transcripts (94.8% of total predicted unique transcripts), corresponding to 20,177 unique descriptions (UD). Differential gene expression analysis resulted in 10,859 transcripts significantly changing among groups when One-way ANOVA (P < 0.05) was applied. The number of differentially expressed (DE) transcripts decreased to 2,958 with FDR-adjusted P < 0.05, being used this set of genes for initial comparisons among groups. Such approach yielded 2,009 DE transcripts (1324 UD) when comparing the diets (FUTURE vs CTRL) within the GS lineage, whereas the same comparison in REF fish only disclosed 323 DE transcripts (207 UD) (Fig. 1). On the other hand, similar results were found when comparing genotypes (GS vs REF) within the CTRL-fed fish (406 DE, 263 UD) or within the animals receiving the FUTURE diet (372 DE, 233 UD). These findings indicate a strong interaction between nutritional programming and genetic background when matching the differential offspring transcriptional response to diet for a given genotype, but not for a given diet and the achieved response with REF and GS fish genotypes. As a validation procedure, six genes covering a wide range of up- and down-regulation between FUTURE and CTRL diets in GS fish were selected, and their fold-change values calculated using real-time PCR were quite consistent (r = 0.998) with those of the RNA-seq analysis (Additional file 2: Supplementary Table 2).

Fig. 1
figure 1

Differentially expressed (DE) transcripts in genetically selected fish (GS) and reference fish (REF) fed either FUTURE or CTRL diets during the challenge phase. Numbers indicate differentially expressed transcripts (FDR-adjusted P < 0.05)

Full size image

Discriminant and Gene Ontology analysis

To better explore the different nutritionally mediated response of the offspring with a different genetic background, a partial least squares-discriminant analysis (PLS-DA) was conducted with the data filtered by One-Way ANOVA in GS fish fed CTRL or FUTURE diets, and REF fish fed CTRL or FUTURE diets (Fig. 2). The two discriminant models were validated significantly for either GS fish (Additional file 3: Supplementary Fig. 1A) or REF fish (Additional file 3: Supplementary Fig. 1B), being all animals correctly classified in each group by hierarchical cluster analysis (Additional file 3: Supplementary Figs. 1D and E). However, the number of discriminant transcripts (VIP ≥ 1) was greater in GS fish than in REF fish (3,648 vs 2,909). Indeed, the explained variance was similar in the two genetic groups of fish (R2X = 93%-99%), whereas the predicted variance was higher in GS fish (Q2 = 74%) than in REF fish (Q2 = 21%) (Fig. 2).

Fig. 2
figure 2

Scores plot of partial least-squares discriminant analysis (PLS-DA) of hepatic transcripts after the challenge phase with CTRL and FUTURE diets in GS (A) and REF (B) fish. RNA-seq data in the analysis were normalized values of differentially expressed transcripts (One-way ANOVA, P < 0.05). The number of discriminant transcripts (VIP ≥ 1; in grey), together with the filtered and enriched functions (GO-BP; in purple) are shown at the bottom right of each plot

Full size image

The over-representation analysis of the list of discriminant transcripts discerned a larger number of enriched functions (Gene ontology—Biological process, GO-BP, unique terms) that was higher in GS fish (546) than in REF fish (242) (Additional file 4: Supplementary Tables 3A and B). For the simplification of the analysis, data were filtered for the enriched functions of upper levels of GO-BP categories, reducing the numbers of terms to 185 in GS fish and 118 in REF fish (Figs. 2A and B). The resulting over-represented functions were clustered in 39 supra-categories (GO-BP ancestors), and the different numbers of DE transcripts within each one are presented in Figs. 3A and C. In both GS and REF fish, the GO-BP supra-categories Localization and Response to stimulus showed the largest number of DE transcripts, although other relevant supra-categories, such as Lipid metabolic process or N compound metabolic process were largely represented. Interestingly, most of these supra-category transcripts (nearly 92%) were down-regulated by FUTURE diet in the offspring of challenged GS fish, whereas the opposite trend occurred in REF fish (nearly 65% of all transcripts were up-regulated). Networks were also performed to show the associations among the different supra-categories of over-represented GO-BP terms, which in turn resulted in an extensive number of links in both GS (Fig. 3B) and REF fish (Fig. 3D).

Fig. 3
figure 3

Bar plot depicting the results of an over-representation test performed over the GO-BP terms of the filtered transcripts for the GS (A) and REF (C) fish. These transcripts were classified in the different GO-BP ancestors presented in the figures. The size of the bars represents the number of transcripts, which are up-regulated (in black) or down-regulated (in grey). Network layout representing the associations between the assigned GO-BP ancestors according to their shared allocated genes in the GS (B) and REF (D) fish. Node size represents the number of transcripts and node colours, the representative name of GO-BP ancestor. Edge width represents the number of shared genes between two supra-categories. * indicates that the supra-category appears in both GS and REF fish. Met., metabolic; Cel., cellular; Reg., regulation; Multicel. org., multicellular organismal

Full size image

Patterns of offspring DNA methylation

Using MBD-seq, ~ 750 million single-end (SE) reads (1 × 75), containing at least one methylated CpG, were obtained, with an average of ~ 31 million reads per sample (Additional file 1: Supplementary Table 1). After trimming and quality filtering, around 6% of all liver reads were discarded, and the remaining reads ranged between 22 million (~ 1.65 Gb) and 40.5 million (~ 3 Gb) in all samples. Up to 84% of these pre-processed reads were mapped against the CSIC gilthead sea bream reference genome (1.6 Gb), which was divided in 25-bp windows (the unit for calculating the level of methylation). Such analysis identified a total of ~ 10 M 25-bp genomic regions that appeared to contain CG dinucleotides susceptible to be altered by methylation. These methylated regions spanned ~ 263 Mb of the total gilthead sea bream genome, and comprised ~ 13.6 M CpG (74.3% of the total 18.3 M CG in the gilthead sea bream genome).

Differential methylation studies detected a total of ~ 553,000 25-bp genomic regions changing (P < 0.05), at least, in one comparison within the four experimental groups. The normalized methylation values (rpkm, reads per kilo base per million mapped reads) of these differentially methylated (DM) regions were used as input in discriminant analysis to assess the effect of nutritional programming and genetic background over DNA-methylation patterns. Unlike RNA-seq results, when both dietary groups were compared within each genetic background, the discriminant separation was not statistically significant (P > 0.05). Nonetheless, when the four groups were analysed together, the discriminant model was statistically validated (Additional file 3: Supplementary Fig. 1C), being 21 out of 24 fish correctly classified by hierarchical clustering analyses (Additional file 3: Supplementary Fig. 1F). The resulting model showed percentage values of explained (R2X) and predicted (Q2) variance that remained above 60% and 38%, respectively (Fig. 4). When VIP threshold (VIP ≥ 1) was applied in the validated PLS-DA model, up to 223,154 25-bp DM genomic regions drove the separation between experimental groups. These DM regions of discriminant value were found to be located in the promoter/coding sequence of over 34,000 DM transcripts.

Fig. 4
figure 4

Scores plot of partial least-squares discriminant analysis (PLS-DA) of methylated regions in GS and REF fish fed the FUTURE (in red) or the CTRL (in black) diet after the challenge phase. MBD-seq data were the normalized values of differentially methylated 25-bp genomic regions (One-way ANOVA, P < 0.05). Graphical representation of the goodness-of-fit of the PLS-DA model is presented at the upper left of the plot

Full size image

Differential gene expression and DNA methylation overlapping

When the previously mentioned discriminant DE and DM transcripts were matched, a total of 1,785 (1,592 UD) and 1,548 (1,397 UD) transcripts were overlapping in GS (Fig. 5A) and REF (Fig. 5B) fish, respectively. Functional network analysis of these overlapping genes of fish receiving FUTURE or CTRL diets displayed 195 (854 transcripts; 751 UD; 4,975 DM regions) and 42 (383 transcripts; 327 UD; 2,856 DM regions) different GO-BP terms in GS (Fig. 5A) and REF (Fig. 5B) fish, respectively (Additional file 4: Supplementary Tables 3C and D). Then, we kept with those DM regions that showed an opposite trend for DNA-methylation and expression (i.e., hyper-methylated regions with down-regulation of the matching transcript or vice versa), and after filtering for transcripts with at least 80% of their DM regions with a negative correlation, we displayed a total of 2,348 and 805 DM regions in GS fish (Fig. 6A) and REF fish (Fig. 6C), respectively. Of note, DM regions corresponding to 30 GO-BP ancestors were located in 264 (246 UD) and 99 (96 UD) DE transcripts of GS (Fig. 6B) and REF fish (Fig. 6D), respectively. Most of these DM regions (nearly 88%) were hyper-methylated in GS fish fed FUTURE diet compared to those receiving the CTRL diet, and subsequently their matching genes were down-regulated. Conversely, the opposite pattern occurred in REF fish, where around 66% of all DM regions were hypo-methylated in animals fed FUTURE diet vs those ingesting the CTRL diet. The distinct effect of diet depending on the genetic background seems to indicate again a clear interaction between genetic lineage and nutritional programming. The wide distribution of DM regions across introns (48%), exons (27.4%) and promoters (24.6%) is shown for GS and REF fish in Additional file 5: Supplementary Fig. 2.

Fig. 5
figure 5

Venn diagrams showing the results of overlapping differentially methylated (DM) transcripts by MBD-seq (One-way ANOVA, P < 0.05, VIP ≥ 1) and differentially expressed (DE) transcripts by RNA-seq (One-way ANOVA, P < 0.05; VIP ≥ 1) following the challenge phase with FUTURE and CTRL diets in GS (A) and REF (B) fish. The numbers of unique gene descriptions corresponding to the overlapping genes are in blue. The numbers of enriched functions (GO-BP terms, in green) and their unique gene descriptions (in red) are also presented

Full size image
Fig. 6
figure 6

Correlation analysis between differentially methylated (DM) genomic regions and their corresponding transcripts (both as log2 fold change in rpkm values between FUTURE and CTRL-fed fish) in GS (A) and REF (C) fish. Significance of correlation is shown. The results of an over-representation test performed over the GO-BP terms of filtered DM regions showing opposite trends for DNA methylation and the expression of the matching gene are also presented for GS (B) and REF (D) fish. The GO-BP ancestors assigned to these selected transcripts are presented. The size of the bars represents the number of DM regions, which are hyper-methylated (in black) or hypo-methylated (in white). * indicates that the supra-category appears in both GS and REF fish

Full size image

Candidate offspring epigenetic markers

Among the 30 selected GO-BP ancestors, we filtered those (11) with over 85% of genes allocated to the enriched GO-BP terms of GS fish (Fig. 5A), and not to the enriched ones of the REF fish (Fig. 5B). For that, the top-five GO-BP ancestors ordered by the number of genes were Cell cycle, Viral process, Protein transport, Lipid metabolic process, and Cellular process (Fig. 7A). These path ancestors allocated 115 transcripts (106 UD) (Additional file 6: Supplementary Table 4), which were further ordered by their expression fold-change when comparing challenged FUTURE and CTRL dietary groups (Fig. 7B). As a result of this, Lipid metabolic process showed the greatest percentage of genes within the first quartile (~ 40%), followed by Cell cycle (~ 25%). The rest of GO-BP ancestors (Protein transport, Cellular process and Viral process) remained below 20%. Regarding the in-depth study of the 23 genes related to Lipid metabolic process, the position of the DM regions and the number of CpGs are presented in Table 2 and Additional file 7: Supplementary Fig. 3. The genomic organization and number of position of DM CpGs was graphically represented for the top-ten genes (cd36, fatty acid translocase; cidea, cell death-inducing DNA fragmentation factor; fasn, fatty acid synthase; g6pd, glucose-6-phosphate dehydrogenase; acsbg2, long-chain-fatty-acid-CoA ligase; pitpna, phosphatidylinositol transfer protein alpha; acsbg2, long-chain-fatty-acid-CoA ligase; lipt1, lipoyltransferase; scd1a; acsl4, acyl-CoA synthetase 4) within the category of Lipid metabolic process (Fig. 8A). Within this top list, those transcripts with a DM region spotted in the promoter presented a higher average value of Log2 fold-change in their expression (~ 2.5) than those with their DM regions in exon (~ 1.6) or intron (~ 1.5) genomic regions (Fig. 8B). When searched for enriched transcription factor binding sites (TFBS), hepatocyte nuclear factor 3-beta (HNF-3) and transforming growth factor-beta 1 (AP-1) appeared, but none of them were located within the DM regions of any of the 23 genes related to Lipid metabolic processes (Fig. 8A, Additional file 7: Supplementary Fig. 3).

Fig. 7
figure 7

Bar chart representing the top five GO-BP supra-categories of GS fish showing a negative correlation for at least the 80% of differentially methylated regions and their associated differentially expressed transcripts (A). Green colour indicates hyper-methylation (FUTURE vs CTRL fed fish), whereas red means hypo-methylation. List of all the genes represented in A are ordered by the value of log2 fold-change in gene expression, together with a bar chart with the percentage of genes falling within the first quartile of the mentioned list (considering the absolute value of the log2 fold-change) (B). In the ordered list, each colour represents a different supra-category, being in blue the genes of the supra-category Lipid metabolic process. The down-regulated genes appear above the line on the left (indicated by a blue arrow) and the up-regulated genes below (green arrow)

Full size image
Table 2 Differential expression (log2 fold-change FUTURE vs CTRL diet, log2 FC (exp)) and location (scaffold and chromosome, chr) of the selected 23 genes (potential epigenetic markers) belonging to the supra-category Lipid metabolic process and their corresponding putative differentially methylated (DM) regions. The genomic area (Feature: promoter, Prom, exon or intron) where these DM regions are positioned and their differential methylation (log2 fold-change FUTURE vs CTRL diet, log2 FC (met)) are also shown
Full size table
Fig. 8
figure 8

Genomic organization of the 10 top genes belonging to the GO-BP Lipid metabolic process with the greatest absolute change in gene expression during the challenge phase in genetically selected fish (A). Red boxes indicate the location of the differentially methylated regions. The number of CpGs within each region is shown below the gene representation and the coloured squares refer to the presence of transcription factor binding sites (green for AP-1 and blue for HNF-3). Bar chart representing the average absolute value of log2 fold-change in gene expression (FUTURE vs CTRL) for the genes shown in A whose differentially methylated regions appear in promoter, exon or intron genomic areas (B)

Full size image

Discussion

Evidences showing a link between maternal and early life nutrition and epigenetics have been shown and reviewed by several authors in farmed fish [56, 57]. Hence, the transcriptional impact of epigenetic modifications has an effect on the organism's phenotype that results in developmental adaptations that permanently change the structure, physiology and metabolism of affected animals during adult life [58]. Thus, it is well recognized that environmental regulation of economically important aquaculture traits (e.g., growth, disease resistance, low oxygen and warmer temperature tolerance, etc.) is mediated, at least in part, at the level of epigenetic regulation, and such environment-induced epigenetic changes appeared to be intergenerationally inherited, though evidences for transgenerational inheritance are still limited [59]. The fine interplay of genetics and epigenetics is thereby hardly underlying, though a recent study highlighted that the global hepatic methylome landscape and the expression level of DNA (de)methylation-related genes were differentially regulated in trout isogenic lines [60]. In the same line, we found herein that nutritional programming differentially affects the hepatic transcriptome and genome-wide DNA methylome of farmed gilthead sea bream depending on genetic background. Certainly, the offspring of GS fish within the PROGENSA® breeding program shared a better performance than those of REF animals during the challenge phase, and differences observed between fish fed with the CTRL or FUTURE-based diet seemed to be lower in GS lineage. This nutritional challenge was in turn associated with a great impact on the number of genes and functions potentially altered by epigenetic mechanisms, especially in GS animals. Such statement about the latter fish was especially evident for the GO-BP term Lipid metabolic process after the application of different filters that ultimately reflected in a consistent manner the offspring epigenetic-mediated changes on hepatic lipogenesis and fatty acid metabolism. This, together with the observation that, during the challenging phase, the number of DE genes and the associated GO-BP terms was larger in GS fish than in REF fish highly supported an improved metabolic plasticity of GS lineage (Figs. 1, 2 and 3), which was first evidenced herein by a wide-liver transcriptome response. Likewise, previous microbiota studies also supported an improved resilience and functional plasticity of GS fish that can become especially relevant in a scenario of climate change. Certainly, using metatranscriptomic and inferred-metagenome approaches, we found that the adherent intestinal microbiota of GS fish may change their function and activity instead of their composition to cope with changes in diet composition. Such statement was first evidenced at harvest with the offspring of F1/F2 crosses farmed in a Mediterranean experimental facility [53, 54], and thereafter this has been confirmed on a temporal basis through the production cycle of F3/F4 crosses at Canary Islands [55].

Concerning DNA methylation, a large body of evidence in a wide range of fish links changes in epigenetic signatures with a different genetic background [59, 61, 62]. This can also apply to gilthead sea bream and accordingly, in the present study, important differences in the methylation patterns were found when comparisons were made between GS and REF fish (Fig. 4). Moreover, when DM regions were matched with DE transcripts, a remarkable number of DM regions corresponding to DE transcripts were found in response to nutritional programming in either GS or REF animals (Fig. 5), which is in agreement with previous studies in the same [24, 63] or other fish species [8, 10, 64]. However, the number of GO-BP functions associated to DE transcripts whose expression could be potentially altered by different DNA methylation patterns was over 4-fold larger in GS than in REF fish (Fig. 5), which might also imply a greater functional plasticity of the GS linage in terms of epigenetic regulatory mechanisms.

Methylated CpGs are usually associated with the silencing of gene expression [65], so we retained those transcripts showing negative correlation between DNA methylation and their expression level in FUTURE vs. CTRL-fed animals within each fish lineage (Fig. 6). Interestingly, most of these DM regions (approximately 88% of all) were hyper-methylated in GS fish receiving the FUTURE diet, whereas in REF animals most of them (approximately 66% of all) were hypo-methylated. These percentages of hyper and hypo-methylation were similar to those calculated for down (almost 92%) or up-regulated (nearly 65%) DE transcripts in the offspring of GS and REF fish during the challenge phase (Fig. 3). Such analogy further supports the relevance of DNA methylation in the regulation of gene expression in fish and in our experimental model in particular [9]. Otherwise, like in humans and other animal models [66, 67], a global DNA hypo-methylation and site-specific hyper-methylation could be associated with aging in marine fish [68, 69]. We may then speculate that changes in the epigenetic clock of REF animals might occur, reflecting perhaps an aging phenotype with an overall impairment of performance that affects among other traits, growth, swimming behaviour and the incidence of skeletal deformities [49,50,51,52]. All this reinforces the role of differentially DNA methylation patterns as key epigenetic markers, and after applying stringent filters focused on the differential response of GS fish, up to 23 genes with the GO-BP term Lipid metabolic process were retained as powerful epigenetic markers of early nutritional programming of hepatic metabolism with FUTURE-based diets. An initial prospect of DM regions revealed their wide distribution across the genome, but it is usually believed that differentially methylated CpG sites in the promoter region exert a stronger regulatory effect upon gene expression [70]. Accordingly, when responsive genes related to lipid metabolism were ordered by fold-change in expression, the top 10 genes showed a higher concentration of differentially methylated CpG sites in their promoter region (Fig. 8B). Nonetheless, it should be born in mind that the relationship of the epigenetic modifications with gene expression seems to be far more complex than initially thought [71], and regulatory DM regions in areas different to promoters have been reported across all vertebrate species [31, 72, 73].

Regarding in-depth the list of the top 10 candidate epigenetic markers, scd1a was among the selected genes. This gene encodes for the stearoyl-CoA desaturase, which has a role in the synthesis of 18C and 16C monounsaturated fatty acids (oleic acid and palmitoleic acid) from stearoyl-CoA and palmitoyl-CoA, respectively [74]. In a previous gilthead sea bream study, a reduced expression of scd1a was concurrent with a higher DNA methylation level in the gilthead sea bream offspring of parents fed an ALA-enriched diet with a limited supply of n-3 LC-PUFA [24]. In particular, that regulation was associated to an increased methylation at a regulatory region in the proximal scd1a promoter, a genomic area rich in hypo-methylated CGI [75]. Conversely, in the present study a differentially methylated cytosine was found in the second intron instead of the promoter region (Fig. 8A). It can be argued that fish and experimental conditions are not the same in this and the previous study, but importantly Perera et al. [24] focused on CGI of the proximal promoter that mostly becomes hypo-methylated, especially in the case of genes with a constitutive high expression level [76]. By contrast, we used herein a wide-genome approach (MBD-seq) that primes the recognition of differentially methylated CpG sites on methylated genomic regions [41]. In this regard, the no coincidence of DM regions in this and the previous [24] study could be mainly attributed to the complementarity of methodologies rather than to the discordance of results.

Besides scd1a, other candidate epigenetic markers are the g6pd gene that encodes for the glucose-6-phosphate dehydrogenase, a major regulatory enzyme involved in the generation of NADPH, which is required by the fasn for catalysing all the reaction steps involved in the conversion of acetyl-CoA and malonyl-CoA to palmitate [77, 78]. Other important lipogenic enzyme is lipt1, a lipoyltransferase that is required for 2-ketoacid dehydrogenase function and mitochondrial fatty acid synthesis [79]. The relative contribution of all these genes, in addition to that of scd1a, on the resulting lipid metabolism phenotype remains elusive, but all of them showed significantly hyper-methylated CpG sites with a concurrent down-regulated expression in the offspring of GS fish fed the FUTURE diet. Similar DNA-methylation and gene expression patterns were found for cidea and pitpna. The former is a key landmark of apoptosis that plays a crucial role in lipid and energy metabolism including lipolysis, lipid oxidation and lipid droplet formation, resulting its low expression in a reduced accumulation of lipid depots [80]. Likewise, the pitpna gene encodes for a lipid-binding protein that catalyzes the transfer of phosphatidylinositol and phosphatidylcholine from the Golgi apparatus to the endoplasmic reticulum [81], and its up-regulated expression in zebrafish is concurrent with signs of hepatic steatosis [82]. Overall, this gene expression and epigenetic down-regulated pattern would drive the hepatic lipid metabolism towards reduced hepatic fatty acid biosynthesis and lipid storage, which will act limiting the lipogenic pathway in response to a low fish oil diet during the challenge phase [24]. The fact that this outcome appeared in GS animals is consistent with a recent study where the expression of several hepatic gene markers of lipid metabolism, including scd1a, were down-regulated in gilthead sea bream differentially selected for improved feed conversion ratio, suggesting that more efficient fish are also likely to present lowered hepatic lipogenesis and fat deposition rates [83].

Acyl-CoA synthetases, encoded by acsl4 and acsbg2 genes, were also identified as robust epigenetic markers of nutritional programming in our model of GS fish. These genes encode for enzymes catalysing the conversion of long-chain fatty acids to their active form acyl-CoA for both synthesis of cellular lipids and degradation via β-oxidation. The Acsl4 enzyme shows preference for arachidonic acid and EPA as substrates [84], whereas the latter has increased ability for oleic acid and linoleic acid [85]. We may then speculate that, in our experimental conditions, with a low fish oil diet, the enhanced DNA methylation and down-regulated expression of acsl4 and acsbg2 genes will serve to limit and preserve the use of unsaturated fatty acids for vital functions of the organism. This would be consistent with the idea of an enhanced mobilization of LC-PUFA from liver to other tissues in response to similar attempts of nutritional programming in the same species [23]. Lastly, among the 10 top selected genes, the only one showing hypo-methylation and up-regulation (FUTURE vs CTRL) was cd36, whose encoded protein is an important fatty acid transporter associated with long-chain fatty acid trans-membrane uptake [86,87,88]. This would favour the uptake of available fatty acids, which would be compatible with the described down-regulation of genes related to lipogenesis and hepatic fat storage. Accordingly, Turkmen et al. [13] allowed that early nutritional programming with linseed oil-rich diets reduced the risk of hepatic steatosis in gilthead sea bream. Overall, all this agrees with the initial statement that the epigenetic regulation of gene expression due to nutritional programming works towards a balanced physiological response of the animals by means of precluding an over-expression of specific genes that might result counterproductive in a changing environment [24]. Likewise, experimental evidences in salmon indicated that this general view not only applies to nutritional programming, but also to DNA methylation dynamics of fish challenged with high temperature and moderate hypoxia in a context of global warm [76]. Whether those induced epigenetic marks (mediating the response to nutritional programming) might be transmitted over generations cannot be assessed in this study and further research about their transgenerational inheritance needs to be addressed for its application in aquaculture practice.

Conclusions

Gene expression profiles and DNA methylation signatures following nutritional programming were clearly dependent on genetic background in our experimental model. Such assumption affected the magnitude, but also the type and direction of change when comparing GS and REF fish. Accordingly, the resulting epigenetic clock of REF fish might depict an older phenotype with a lower DNA methylation state of DE transcripts during the challenging phase with the FUTURE-based diet formulation. Therefore, attempts to search and validate robust epigenetic markers will be specific of each lineage, and focusing on GS fish we successively applied different filters summarized in Fig. 9, which allows us to reduce the search to 115 candidate epigenetic markers. However, genes with the GO-BP term Lipid metabolism were markedly reordered after filtering by fold-changes in expression. This rendered a final list of top 10 epigenetic markers ordered by the magnitude of response (cd36 > cidea > fasn > g6pd > acbsg2 > pitpna > acsbg2 > lipt1 > scd1a > acsl4), which reinforces the key role of maintaining regulated hepatic lipogenesis and fatty acid metabolism when precluding the over-expression of specific genes that might result counterproductive in a changing environment.

Fig. 9
figure 9

Chart representing the filters and criteria applied for matching wide-transcriptomics (RNA-seq) and wide-genomic DNA methylation (MBD-seq) approaches for recording robust epigenetic markers of special relevance to evaluate the success of nutritional programming in genetically selected (GS) gilthead sea bream within the PROGENSA® program. For the 10 top genes identified as candidate markers, blue colour indicates down-regulation in GS fish challenged with the FUTURE diet, whereas green refers to up-regulation

Full size image

Methods

Ethics statement

All procedures of animal rearing and sampling were carried out according to European animal directives (2010/63/EU) and Spanish laws (Royal Decree RD53/2013) for the protection of animals used in scientific experiments. The Bioethical Committee of the University of Las Palmas de Gran Canaria approved all the protocols used in the present study (approval no. OEBA_ULPGC_26/2019).

Broodstock crosses

A population of 6,122 adult fish from the Canary Islands at the 3rd generation of National Breeding Program (PROGENSA®) were evaluated for growth. The estimated breeding values (EBV, expressed as g of whole body) ranged between -159.14 for reference fish (REF) and + 223.18 for the selected fish (GS) with an average value of 8.59 and a standard deviation value of 52.84. A subset of 196 fish (98 fish per broodstock) was then selected as breeders with values for the EBV varying from -25.95 in the group of REF fish to + 39.68 in the group of GS fish, comprising almost the 47% of the evaluated population.

Diets and experimental design

Two experimental, with a varying granule size and composition, and the broodstock diets were employed in this study. The control diet (CTRL) contained fish meal (15%) as the main protein source. The alternative diet (FUTURE) was half-reduced in fish meal (7.5%). Fish oil was added at a relatively high level (5.7–7.6%) in CTRL diet, whereas the FUTURE diet was completely devoid of fish oil. Ingredients and composition are presented in Tables 3 and 4. Gilthead sea bream brood fish were used in this study, belonging to fish genetically selected for high growth (GS) or to Reference fish (REF) within the PROGENSA® selection program. Briefly, both groups, GS and REF, were kept in 1,000 L tanks in a flow-through system with filtered seawater, strong aeration, natural photoperiod and temperature conditions in Canary Islands latitude (27º 59’ N; 15º 22’ W) in the experimental facilities of IU-ECOAQUA (University of Las Palmas de Gran Canaria, Spain). All broodstock animals were fed a low fish oil diet for several months prior to the spawning season (stimulus phase; Fig. 10). Ingredients and composition of the broodstock diet are given in Additional file 8: Supplementary Table 5. Eggs were collected at spawning from the two broodstocks and the offsprings were reared under standard conditions for approximately 5 months (intermediate phase). For further details on the broodstock diet and the spawning quality please see Ferosekhan et al. [89]. Then, 1,600 fish (approximately 12 g weight) from the two broodstock groups (800 fish per genotype) were distributed in 24 500 L tanks and fed (6 days per week) either the FUTURE or the CTRL diet twice a day (08:00 and 14:00 h) for approximately 6 months (challenge phase; Fig. 10). At the end of the trial, 24 h-fasted juvenile fish (6 per diet) were anaesthetized using an overdose of natural clove oil. Livers were rapidly excised, frozen in liquid nitrogen and stored at − 80 °C until RNA and DNA extraction for analyses of gene expression and DNA methylation, respectively.

Table 3 Ingredients and proximate composition of the experimental diets (control, CTRL, and FUTURE) with different pellet sizes (1.8 and 3 mm)
Full size table
Table 4 Fatty acid composition (% of Total Fatty Acids—TFA) of the experimental diets (control, CTRL, and FUTURE)
Full size table
Fig. 10
figure 10

Diagram of the experimental design

Full size image

DNA and RNA Extraction

Liver DNA was extracted using the Quick-DNA™ Mini-prep Plus Kit (Zymo Research, Irvine, CA, USA) following the manufacturer’s instructions. The quantity and quality of DNA were assessed by a NanoDrop 2000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and DNA integrity was assessed in 1% agarose gels. Total RNA (70–100 μg) from liver was extracted with the MagMAX™-96 Total RNA Isolation Kit (Applied Biosystems, Foster City, CA, United States). The RNA concentration and purity was determined using a NanoDrop 2000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Quality and integrity of the isolated RNA were checked on an Agilent Bioanalyzer 2100 total RNA Nano series II chip (Agilent, Amstelveen, Netherlands), yielding RNA integrity numbers (RINs) between 8 and 10. Samples were stored at -80 ºC until DNA and RNA sequencing.

DNA and RNA Illumina sequencing

For the methyl-CpG-binding domain sequencing (MBD-seq) analysis, 300 ng of DNA were fragmented to 200–550 bp using the methylation-insensitive restriction enzyme MseI (New England Biolabs, United States), which recognizes genomic T↓TAA sites, typically found outside of CGIs. The enzyme action and inactivation temperatures, as well as its actuation times and concentration were fixed according to manufacturer’s indications. Products were obtained using AMPure beads and checked and quantified with Picogreen (Invitrogen, Carlsbad, United Stated). Fragmented DNA was then submitted to methylation enrichment using the MethylCollector™ Ultra kit (Active Motif, Carlsbad, CA, United States), following the instructions. Briefly, methylated DNA was captured from 75 ng of fragmented DNA via binding to the methyl-CpG binding domain of the MBD2 protein. Illumina MBD-seq libraries were prepared from 15 ng of methylated DNA fragments using the NEBNext® Ultra™ II DNA Library Prep Kit (Illumina Inc. San Diego, CA, USA) according to the manufacturer’s instructions. All libraries were sequenced on an Illumina NovaSeq 6000 sequencer as a 1 × 75 nucleotides SE read format, according to the manufacturer’s protocol. Raw sequenced data were deposited in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the Bioproject accession number PRJNA915228 (BioSample accession numbers: SAMN32381315-SAMN32381338).

Illumina RNA-seq libraries were prepared from 500 ng total RNA using the Illumina TruSeq™ Stranded mRNA LT Sample Prep Kit (Illumina Inc. San Diego, CA, USA) according to the manufacturer’s instructions. All RNA-seq libraries were sequenced on an Illumina NovaSeq 6000 sequencer as 2 × 150 nucleotides paired-end (PE) read format according to the manufacturer’s protocol. Raw sequenced data were deposited in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the Bioproject accession number PRJNA915228 (BioSample accession numbers: SAMN32381339-SAMN32381362).

Bioinformatics analyses

After sequencing, the quality of the RNA-seq and MBD-seq resulting raw reads was evaluated with FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ accessed on 16 April 2020). In the case of RNA-seq data, libraries were filtered with Trimmomatic (RRID:SCR_011848) [90] eliminating those with quality < 18, length < 200 bp, and > 5% of Ns in the sequence. Cleaned reads were mapped against gilthead sea bream reference genome [91] (available at https://seabreamdb.nutrigroup-iats.org/), using STAR (RRID:SCR_004463) [92]. Unique transcript hit counts were calculated by using featureCounts from the Subread package [93].

In the case of MBD-seq data, pre-processing was performed with Prinseq (RRID:SCR_005454) [94], eliminating those with quality < 26, length < 60 bp, and > 5% of Ns in the sequence. Before mapping, the repeat regions of the gilthead sea bream reference were masked using the BSgenome R package. Then, high-quality reads were aligned to this masked genome using the Bowtie2 software (RRID:SCR_016368) [95]. After MBD-seq mapping, the methylated reads were located in their corresponding genomic region using five training sets created ad-hoc for this analysis. The exon and intron training sets contained the coordinates of these elements along the genome. The promoter region of a gene was considered 5,000 bp before its start codon. Finally, two training sets were established for CGI, depending on if these islands were found in promoters or intergenic regions. Predictions of CGI were done with NewCpGReport tool from the EMBOSS suite [96]. Search parameters for CGI were: length ≥ 200, C + G content ≥ 50%, ratio of observed/expected CpGs ≥ 0.60 and window size = 100.

Statistical analysis

Data of fish performance were analysed using the analysis of variance (ANOVA) at a significance level of 0.05. Two-way ANOVA was applied to these results to determine the combined effects of genotype (GS or REF), diet (FUTURE or CTRL) and their interaction, using the program IBM SPSS version 20 for Windows (IBM SPSS Inc., Armonk, NY, USA).

DE transcripts in RNA-seq data were retrieved using DESeq2 at two significance thresholds (P < 0.05 and FDR < 0.05) [97]. DM regions in MBD-seq data were obtained using the MEDIPS R package (P < 0.05) [98] over regions of 25-bp size, containing at least one dinucleotide CG. To study the separation of the experimental groups, we performed several partial least-squares discriminant analysis (PLS-DA) using EZinfo v3.0 (Umetrics, Umeå, Sweden). DE transcripts and DM regions with a P < 0.05 were introduced in the analyses. The fitness and predictability of these models were validated by a 500 random permutation test (pR2Y < 0.05; pQ2 < 0.05) using the ropls R package [99]. The discriminant ability of each marker was ranked after the creation of the models according to its Variable Importance in the Projection (VIP), and useful markers were detected under a VIP ≥ 1 [100]. Over-representation tests of GO-BP terms and TFBS were implemented in the goseq R package [101] and statistical significance was accepted at FDR < 0.05. GO-BP levels and supra-categories were retrieved using GOATools [102]. All the networks in this work were performed with Cytoscape (RRID:SCR_003032) v2.8 [103]. Pearson correlation coefficients between the expression and methylation Log2FC were performed using the SigmaPlot v14.5 (Systat Software Inc.) software. Gene structure representations were obtained using the genemodel R package (https://github.com/greymonroe/genemodel) and the IBS software [104].

Availability of data and materials

Most data generated or analysed during this study are included in this article and its supplementary information files. Raw sequenced data were deposited in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the Bioproject accession numbers PRJNA915228 (BioSample accession numbers: SAMN32381315-SAMN32381338) and PRJNA915228 (BioSample accession numbers: SAMN32381339-SAMN32381362).

Abbreviations

ALA:

α-linolenic acid

ANOVA:

Analysis of variance

AP-1:

Transforming growth factor-beta 1

BW:

Body weight

CGI:

CpG island

CTRL:

Control diet

DE:

Differentially expressed

DGI:

Daily growth index

DHA:

Docosahexaenoic acid

DM:

Differentially methylated

EPA:

Eicosapentaenoic acid

EBV:

Estimated breeding values

FCR:

Feed conversion ratio

FUTURE:

FUTURE diet

GS:

Genetically selected

GO-BP:

Gene ontology Biological process

HNF-3:

Hepatocyte nuclear factor 3-beta

K:

Condition factor

LC-PUFA:

Long-chain polyunsaturated fatty acids

MBD-seq:

Methyl-binding domain sequencing

NCBI:

National Center for Biotechnology Information

PE:

Paired-end

PLS-DA:

Partial least-squares discriminant analysis

PUFA:

Polyunsaturated fatty acids

REF:

Reference

RINs:

RNA integrity numbers

RNA-seq:

Massive sequencing of mRNA

SE:

Single-end

SGR:

Specific growth rate

SL:

Standard length

SRA:

Sequence Read Archive

TFA:

Total fatty acids

TFBS:

Transcription factor binding sites

UD:

Unique description

VIP:

Variable importance in the projection

WGBS:

Whole genome bisulphite sequencing

References

  1. Lucas A. Role of nutritional programming in determining adult morbidity. Arch Dis Child. 1994;71:288–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. McMullen S, Langley-Evans SC, Gambling L, Lang C, Swali A, McArdle HJ. A common cause for a common phenotype: The gatekeeper hypothesis in fetal programming. Med Hypotheses. 2012;78:88–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Michonska I, Łuszczki E, Zielinska M, Oleksy Ł, Stolarczyk A, Katarzyna D. Nutritional programming : History, hypotheses, and the role of prenatal factors in the prevention of metabolic diseases– A narrative review. Nutrients. 2022;14:4422.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Kenyon PR, Blair HT. Foetal programming in sheep – Effects on production. Small Rumin Res. 2014;118:16–30.

    Article  Google Scholar 

  5. Nascimento KB, Castilho M, Andr J, Meneses M, Prezotto LD, Haddad L, et al. Effects of maternal protein supplementation at mid-gestation of cows on intake, digestibility, and feeding behavior of the offspring. Animals. 2022;12:2865.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Panserat S, Marandel L, Geurden I, Veron V, Dias K, Plagnes-Juan E, et al. Muscle catabolic capacities and global hepatic epigenome are modified in juvenile rainbow trout fed different vitamin levels at first feeding. Aquaculture. 2017;468:515–23.

    Article  CAS  Google Scholar 

  7. Skjærven KH, Jakt LM, Fernandes JMO, Dahl JA, Adam AC, Klughammer J, et al. Parental micronutrient deficiency distorts liver DNA methylation and expression of lipid genes associated with a fatty-liver-like phenotype in offspring. Sci Rep. 2018;8:1–16.

    Article  Google Scholar 

  8. Zhu QS, Wang J, He S, Liang XF, Xie S, Xiao QQ. Early leucine programming on protein utilization and mTOR signaling by DNA methylation in zebrafish (Danio rerio). Nutr Metab. 2020;17:1–13.

    Article  CAS  Google Scholar 

  9. Hou Z, Fuiman LA. Nutritional programming in fishes: insights from mammalian studies. Rev Fish Biol Fish. 2020;30:67–92.

    Article  Google Scholar 

  10. Callet T, Li H, Coste P, Glise S, Heraud C, Maunas P, et al. Modulation of energy metabolism and epigenetic landscape in rainbow trout fry by a parental low protein/high carbohydrate diet. Biology. 2021;10:585.

  11. Izquierdo MS, Turkmen S, Montero D, Zamorano MJ, Afonso JM, Karalazos V, et al. Nutritional programming through broodstock diets to improve utilization of very low fishmeal and fish oil diets in gilthead sea bream. Aquaculture. 2015;449:18–26.

    Article  CAS  Google Scholar 

  12. Turkmen S, Hernández-Cruz CM, Zamorano MJ, Fernández-Palacios H, Montero D, Afonso JM, et al. Long-chain PUFA profiles in parental diets induce long-term effects on growth, fatty acid profiles, expression of fatty acid desaturase 2 and selected immune system-related genes in the offspring of gilthead seabream. Br J Nutr. 2019;122:25–38.

    Article  CAS  PubMed  Google Scholar 

  13. Turkmen S, Zamorano MJ, Fernández-Palacios H, Hernández-Cruz CM, Montero D, Robaina L, et al. Parental nutritional programming and a reminder during juvenile stage affect growth, lipid metabolism and utilisation in later developmental stages of a marine teleost, the gilthead sea bream (Sparus aurata). Br J Nutr. 2017;118:500–12.

    Article  CAS  PubMed  Google Scholar 

  14. Tocher DR, Ghioni C. Fatty acid metabolism in marine fish: Low activity of fatty acyl ∆5 desaturation in gilthead sea bream (Sparus aurata) Cells. Lipids. 1999;34:433–40.

    Article  CAS  PubMed  Google Scholar 

  15. Xie D, Chen C, Dong Y, You C, Wang S, Monroig Ó, et al. Regulation of long-chain polyunsaturated fatty acid biosynthesis in teleost fish. Prog Lipid Res. 2021;82:101095.

    Article  CAS  PubMed  Google Scholar 

  16. Li Y, Monroig O, Zhang L, Wang S, Zheng X, Dick JR, et al. Vertebrate fatty acyl desaturase with Δ4 activity. PNAS. 2010;107:16840–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Morais S, Castanheira F, Martinez-rubio L, Conceição LEC, Tocher DR. Long chain polyunsaturated fatty acid synthesis in a marine vertebrate : Ontogenetic and nutritional regulation of a fatty acyl desaturase with Δ4 activity. Biochim Biophys Acta. 2012;1821:660–71.

    Article  CAS  PubMed  Google Scholar 

  18. Fonseca-madrigal J, Navarro JC, Hontoria F, Tocher DR, Martínez-palacios CA, Monroig Ó. Diversification of substrate specificities in teleostei Fads2: characterization of Δ4 and Δ6Δ5 desaturases of Chirostoma estor. J Lipid Res. 2014;55:1408–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kuah M, Jaya-ram A, Shu-chien AC. The capacity for long-chain polyunsaturated fatty acid synthesis in a carnivorous vertebrate : Functional characterisation and nutritional regulation of a Fads2 fatty acyl desaturase with Δ4 activity and an Elovl5 elongase in striped snakehead (Channa striata). Biochim Biophys Acta. 2015;1851:248–60.

    Article  CAS  PubMed  Google Scholar 

  20. Geay F, Wenon FGD, Mellery J, Tinti E, Mandiki SNM, Tocher DR, et al. Dietary linseed oil reduces growth while differentially impacting LC - PUFA synthesis and accretion into tissues in Eurasian perch (Perca fluviatilis). Lipids. 2015;50:1219–32.

    Article  CAS  PubMed  Google Scholar 

  21. Izquierdo MS, Robaina L, Juárez-Carrillo E, Oliva V, Hernández-Cruz CM, Afonso JM. Regulation of growth, fatty acid composition and delta 6 desaturase expression by dietary lipids in gilthead seabream larvae (Sparus aurata). Fish Physiol Biochem. 2008;34:117–27.

    Article  CAS  PubMed  Google Scholar 

  22. Vagner M, Santigosa E. Characterization and modulation of gene expression and enzymatic activity of delta-6 desaturase in teleosts: A review. Aquaculture. 2011;315:131–43.

    Article  CAS  Google Scholar 

  23. Xu H, Ferosekhan S, Turkmen S, Afonso JM, Zamorano MJ, Izquierdo M. High broodstock fads2 expression combined with nutritional programing through broodstock diet improves the use of low fishmeal and low fish oil diets in gilthead seabream (Sparus aurata) progeny. Aquaculture. 2021;535:736321.

    Article  CAS  Google Scholar 

  24. Perera E, Turkmen S, Simó-Mirabet P, Zamorano MJ, Xu H, Naya-Català F, et al. Stearoyl-CoA desaturase (scd1a) is epigenetically regulated by broodstock nutrition in gilthead sea bream (Sparus aurata). Epigenetics. 2020;15:536–53.

    Article  PubMed  Google Scholar 

  25. Ayisi CL, Zhao JL. Fatty acid composition, lipogenic enzyme activities and mRNA expression of genes involved in the lipid metabolism of Nile tilapia fed with palm oil. Turkish J Fish Aquat Sci. 2017;17:405–15.

    Google Scholar 

  26. Christensen KA, Luyer L, Chan MTT, Rondeau EB, Koop BF, Bernatchez L, et al. Assessing the effects of genotype-by-environment interaction on epigenetic, transcriptomic, and phenotypic response in a Pacific salmon. G3. 2021;11:jkab021.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Ragsdale A, Recalde OO, Dutoit L, Besson AA, Chia JHZ, King T, et al. Paternal hypoxia exposure primes offspring for increased hypoxia resistance. BMC Biol. 2022;20:185.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Lin YL, Zhu ZX, Ai CH, Xlong YY, Liu TD, Lin HR, et al. Transcriptome and DNA methylation responses in the liver of yellowfin seabream under starvation stress. Mar Biotechnol. 2022;25:150–60.

  29. Meissner A, Gnirke A, Bell GW, Ramsahoye B, Lander ES, Jaenisch R. Reduced representation bisulfite sequencing for comparative high-resolution DNA methylation analysis. Nuc Acid Res. 2005;33:5868–77.

    Article  CAS  Google Scholar 

  30. Gu H, Bock C, Mikkelsen TS, Jäger N, Smith ZD, Tomazou E, et al. Genome-scale DNA methylation mapping of clinical samples at single-nucleotide resolution. Nat Methods. 2010;7:133–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Anastasiadi D, Codina AE, Piferrer F. Consistent inverse correlation between DNA methylation of the first intron and gene expression across tissues and species. Epigenet Chromatin. 2018;11:37.

    Article  Google Scholar 

  32. Anastasiadi D, Piferrer F. Epimutations in developmental genes underlie the onset of domestication in farmed European sea bass. Mol Biol Evol. 2019;36:2252–64.

  33. Mukiibi R, Peñaloza C, Gutiérrez A, Yáñez JM, Houston RD, Robledo D. The impact of Piscirickettsia salmonis infection on genome-wide DNA methylation profile in Atlantic salmon. Genomics. 2022;114:110503.

  34. Salem M, Al-Tobasei R, Kenney B. Integrated analyses of DNA methylation and gene expression of rainbow trout muscle atrophy conditions. Genes. 2022;13:1151.

  35. Podgorniak T, Dhanasiri A, Chen X, Ren X, Kuan P, Fernandes J. Early fish domestication affects methylation of key genes involved in the rapid onset of the farmed phenotype. Epigenet. 2022;17:1281–98.

    Article  Google Scholar 

  36. Harris RA, Wang T, Coarfa C, Nagarajan RP, Hong C, Downey SL, et al. Comparison of sequencing-based methods to profile DNA methylation and identification of monoallelic epigenetic modifications. Nat Biotechnol. 2010;28:1097–105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Booth MJ, Branco MR, Ficz G, Oxley D, Krueger F, Reik W, et al. Quantitative sequencing of 5-methylcytosine and 5-hydroxymethylcytosine at single-base resolution. Science. 2012;336:934–7.

    Article  CAS  PubMed  Google Scholar 

  38. Yu M, Hon GC, Szulwach KE, Song CX, Zhang L, Kim A, et al. Base-resolution analysis of 5-hydroxymethylcytosine in the mammalian genome. Cell. 2012;149:1368–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Subhash S, Andersson PO, Kosalai ST, Kanduri C, Kanduri M. Global DNA methylation profiling reveals new insights into epigenetically deregulated protein coding and long noncoding RNAs in CLL. Clin Epigenetics. 2016;8:1–15.

    Article  Google Scholar 

  40. Chan RF, Shabalin AA, Xie LY, Adkins DE, Zhao M, Turecki G, et al. Enrichment methods provide a feasible approach to comprehensive and adequately powered investigations of the brain methylome. Nucleic Acids Res. 2017;45:e97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Aberg KA, Chan RF, van den Oord EJCG. MBD-seq - realities of a misunderstood method for high-quality methylome-wide association studies. Epigenetics. 2020;15:431–8.

    Article  PubMed  Google Scholar 

  42. Aberg KA, Chan RF, Xie L, Shabalin AA, van den Oord EJCG. Methyl-CpG-Binding domain sequencing: MBD-seq. In: Tost J, editor. DNA Methylation Protocols. New York: Springer Nature; 2018. p. 171–90.

  43. Neary JL, Perez SM, Peterson K, Lodge DJ, Carless MA. Comparative analysis of MBD-seq and MeDIP-seq and estimation of gene expression changes in a rodent model of schizophrenia. Genomics. 2017;109:204–13.

  44. Namous H, Peñagaricano F, Del Corvo M, Capra E, Thomas DL, Stella A, et al. Integrative analysis of methylomic and transcriptomic data in fetal sheep muscle tissues in response to maternal diet during pregnancy. BMC Genomics. 2018;19:123.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Mao Z, Li T, Zhao H, Wang X, Kang Y, Kang Y. Methylome and transcriptome profiling revealed epigenetic silencing of LPCAT1 and PCYT1A associated with lipidome alterations in polycystic ovary syndrome. J Cell Physiol. 2021;236:6362–75.

    Article  CAS  PubMed  Google Scholar 

  46. Kim DY, Kim JM. Multi-omics integration strategies for animal epigenetic studies - A review. Anim Biosci. 2021;34:1271–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Woods LC III, Li Y, Ding Y, Liu J, Reading BJ, Fuller SA, et al. DNA methylation profiles correlated to striped bass sperm fertility. BMC Genomics. 2018;19:244.

    Article  PubMed  Google Scholar 

  48. Somerville V, Schwaiger M, Hirsch PE, Walser J-C, Bussmann K, Weyrich A, et al. DNA methylation patterns in the round goby hypothalamus support an on-the-spot decision scenario for territorial behavior. Genes. 2019;10:219.

  49. Perera E, Simó-mirabet P, Suk H, Rosell-moll E, Naya-catalá F, De V, et al. Selection for growth is associated in gilthead sea bream (Sparus aurata) with diet flexibility, changes in growth patterns and higher intestine plasticity. Aquaculture. 2019;507:349–60.

  50. León-Bernabeu S, Shin HS, Lorenzo-Felipe A, García-Pérez C, Berbel C, Elalfy IS, et al. Genetic parameter estimations of new traits of morphological quality on gilthead seabream (Sparus aurata) by using IMAFISH_ML software. Aquac Reports. 2021;21:100883.

  51. Perera E, Rosell-Moll E, Martos-Sitcha JA, Naya-Català F, Simó-Mirabet P, Calduch-Giner J, et al. Physiological trade-offs associated with fasting weight loss, resistance to exercise and behavioral traits in farmed gilthead sea bream (Sparus aurata) selected by growth. Aquac Reports. 2021;20:100645.

  52. Perera E, Rosell-Moll E, Naya-Català F, Simó-Mirabet P, Calduch-Giner J, Pérez-Sánchez J. Effects of genetics and early-life mild hypoxia on size variation in farmed gilthead sea bream (Sparus aurata). Fish Physiol Biochem. 2021;47:121–33.

  53. Piazzon MC, Naya-Català F, Perera E, Palenzuela O, Sitjà-Bobadilla A, Pérez-Sánchez J. Genetic selection for growth drives differences in intestinal microbiota composition and parasite disease resistance in gilthead sea bream. Microbiome. 2020;8:168.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Naya-Català F, Piazzon MC, Calduch-Giner JA, Sitjà-Bobadilla A, Pérez-Sánchez J. Diet and host genetics drive the bacterial and fungal intestinal metatranscriptome of gilthead sea bream. Front Microbiol. 2022;13:883738.

  55. Naya-Català F, Piazzon MC, Torrecillas S, Toxqui-Rodríguez S, Calduch-Giner JA, Fontanillas R, et al. Genetics and nutrition drive the gut microbiota succession and host-transcriptome interactions through the gilthead sea bream (Sparus aurata) Production Cycle. Biology. 2022;11:1744.

  56. Granada L, Lemos MFL, Cabral HN, Bossier P, Novais SC. Epigenetics in aquaculture – the last frontier. Rev Aquac. 2018;10:994–1013.

    Article  Google Scholar 

  57. Best C, Ikert H, Kostyniuk DJ, Craig PM, Navarro-martin L, Marandel L, et al. Epigenetics in teleost fish : From molecular mechanisms to physiological phenotypes. Comp Biochem Physiol Part B. 2018;224:210–44.

  58. Gotoh T. Potential of the application of epigenetics in animal production. Anim Prod Sci. 2015;55:145–58.

    Article  CAS  Google Scholar 

  59. Liu Z, Zhou T, Gao D. Genetic and epigenetic regulation of growth, reproduction, disease resistance and stress responses in aquaculture. Front Genet. 2022;13:994471.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Marandel L, Heraud C, Véron V, Laithier J, Marchand M, Quillet E, et al. A plant-based diet differentially affects the global hepatic methylome in rainbow trout depending on genetic background. Epigenetics. 2022;17:1726–37.

    Article  Google Scholar 

  61. Ding YX, He F, Wen HS, Li JF, Qian K, Chi ML, et al. Polymorphism in exons CpG rich regions of the cyp17-II gene affecting its mRNA expression and reproductive endocrine levels in female Japanese flounder (Paralichthys olivaceus). Gen Comp Endocrinol. 2012;179:107–14.

  62. Si Y, He F, Wen H, Li J, Zhao J, Ren Y, et al. Genetic polymorphisms and DNA methylation in exon 1 CpG-rich regions of PACAP gene and its effect on mRNA expression and growth traits in half smooth tongue sole (Cynoglossus semilaevis). Fish Physiol Biochem. 2016;42:407–21.

  63. Perera E, Yúfera M. Effects of soybean meal on digestive enzymes activity, expression of inflammation-related genes, and chromatin modifications in marine fish (Sparus aurata L.) larvae. Fish Physiol Biochem. 2017;43:563–78.

  64. Xu H, Dong X, Ai Q, Mai K, Xu W, Zhang Y, et al. Regulation of Tissue LC-PUFA Contents, Δ6 Fatty Acyl Desaturase (FADS2) Gene Expression and the Methylation of the Putative FADS2 Gene Promoter by Different Dietary Fatty Acid Profiles in Japanese Seabass (Lateolabrax japonicus). PLoS ONE. 2014;9:e87726.

  65. Jones PA. Functions of DNA methylation : islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13:484–92.

    Article  CAS  PubMed  Google Scholar 

  66. Gentilini D, Mari D, Castaldi D, Remondini D, Ogliari G, Ostan R, et al. Role of epigenetics in human aging and longevity : genome-wide DNA methylation profile in centenarians and centenarians’ offspring. Age. 2013;35:1961–73.

    Article  CAS  PubMed  Google Scholar 

  67. Jung M, Pfeifer GP. Aging and DNA methylation. BMC Biol. 2015;13:7.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Anastasiadi D, Piferrer F. A clockwork fish : Age prediction using DNA methylation-based biomarkers in the European seabass. Mol Ecol Resour. 2020;20:387–97.

    Article  CAS  PubMed  Google Scholar 

  69. Piferrer F, Anastasiadi D. Age estimation in fishes using epigenetic clocks: Applications to fisheries management and conservation biology. Front Mar Sci. 2023;10:1062151.

  70. Moore LD, Le T, Fan G. DNA Methylation and Its Basic Function. Neuropsychopharmacology. 2013;38:23–38.

    Article  CAS  PubMed  Google Scholar 

  71. Yu P, Xiao S, Xin X, Song C, Huang W, Mcdee D, et al. Spatiotemporal clustering of the epigenome reveals rules of dynamic gene regulation. Genome Res. 2013;23:352–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Vanderkraats ND, Hiken JF, Decker KF, Edwards JR. Discovering high-resolution patterns of differential DNA methylation that correlate with gene expression changes. Nuc Acid Res. 2013;41:6816–27.

    Article  CAS  Google Scholar 

  73. Schlosberg CE, Vanderkraats ND, Edwards JR. Modeling complex patterns of differential DNA methylation that associate with gene expression changes. Nuc Acid Res. 2017;45:5100–11.

    Article  CAS  Google Scholar 

  74. Hsieh SL, Chang HT, Wu CH, Kuo CM. Cloning, tissue distribution and hormonal regulation of stearoyl-CoA desaturase in tilapia, Oreochromis mossambicus*. Aquaculture. 2004;230:527–46.

  75. Long HK, King HW, Patient RK, Odom DT, Klose J. Protection of CpG islands from DNA methylation is DNA-encoded and evolutionarily conserved. Nuc Acid Res. 2016;44:6693–706.

    Article  CAS  Google Scholar 

  76. Beemelmanns A, Ribas L, Anastasiadi D, Moraleda-prados J, Zanuzzo FS, Rise ML, et al. DNA methylation dynamics in Atlantic salmon (Salmo salar) challenged with high temperature and moderate hypoxia. Front Mar Sci. 2021;7:604878.

  77. Chen Q-L, Luo Z, Liu X, Song Y-F, Liu C-X, Zheng J-L, et al. Effects of waterborne chronic copper exposure on hepatic lipid metabolism and metal-element composition in Synechogobius hasta. Arch Environ Contam Toxicol. 2013;64:301–15.

  78. Tian J, Wen H, Zeng L, Jiang M, Wu F, Liu W, et al. Changes in the activities and mRNA expression levels of lipoprotein lipase (LPL), hormone-sensitive lipase (HSL) and fatty acid synthetase (FAS) of Nile tilapia (Oreochromis niloticus) during fasting and re-feeding. Aquaculture. 2013;400–401:29–35.

  79. Nowinski SM, Solmonson A, Rusin SF, Maschek JA, Bensard CL, Fogarty S, et al. Mitochondrial fatty acid synthesis coordinates oxidative metabolism in mammalian mitochondria. eLife. 2020;9:e58041.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Chen F-J, Yin Y, Tin Chua B, Li P. CIDE family proteins control lipid homeostasis and the development of metabolic diseases. Traffic. 2020;21:94–105.

    Article  CAS  PubMed  Google Scholar 

  81. Garner K, Hunt AN, Koster G, Somerharju P, Groves E, Li M, et al. Phosphatidylinositol transfer protein, cytoplasmic 1 (PITPNC1) binds and transfers phosphatidic Acid. J Biol Chem. 2012;287:32263–76.

  82. Fai Tse WK, Woei Li J, Kwan Tse AC, Chan TF, Hin Ho JC, Sun Wu RS, et al. Fatty liver disease induced by perfluorooctane sulfonate: Novel insight from transcriptome analysis. Chemosphere. 2016;159:166–77.

    Article  PubMed  Google Scholar 

  83. Calduch-Giner J, Rosell-Moll E, Besson M, Vergnet A, Bruant JS, Clota F, et al. Changes in transcriptomic and behavioural traits in activity and ventilation rates associated with divergent individual feed efficiency in gilthead sea bream (Sparus aurata). Aquac Reports. 2023;29:101476.

  84. Golej DL, Askari B, Kramer F, Barnhart S, Vivekanandan-Giri A, Pennathur S, et al. Long-chain acyl-CoA synthetase 4 modulates prostaglandin E2 release from human arterial smooth muscle cells. J Lipid Res. 2011;52:782–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Pei Z, Jia Z, Watkins PA. The second member of the human and murine “bubblegum” family is a testis- and brainstem-specific. J Biol Chem. 2006;281:6632–41.

  86. Torstensen BE, Nanton DA, Olsvik PA, Sundvold H, Stubhaug I. Gene expression of fatty acid-binding proteins, fatty acid transport proteins (cd36 and FATP) and β-oxidation-related genes in Atlantic salmon (Salmo salar L.) fed fish oil or vegetable oil. Aquac Nutr. 2009;15:440–51.

  87. Zhou J, Stubhaug I, Torstensen BE. Trans-Membrane Uptake and Intracellular Metabolism of Fatty Acids in Atlantic Salmon (Salmo salar L.) Hepatocytes. Lipids. 2010;45:301–11.

  88. Fink IR, Benard EL, Hermsen T, Meijer AH, Forlenza M, Wiegertjes GF. Molecular and functional characterization of the scavenger receptor CD36 in zebrafish and common carp. Mol Immunol. 2015;63:381–93.

    Article  CAS  PubMed  Google Scholar 

  89. Ferosekhan S, Sarih S, Afonso JM, Zamorano MJ, Fontanillas R, Izquierdo M, et al. Selection for high growth improves reproductive performance of gilthead seabream Sparus aurata under mass spawning conditions, regardless of the dietary lipid source. Anim Rep Sci. 2022;241:106989.

    Article  CAS  Google Scholar 

  90. Bolger AM, Lohse M, Usadel B. Genome analysis Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Pérez-Sánchez J, Naya-Català F, Soriano B, Piazzon MC, Hafez A, Gabaldón T, et al. Genome sequencing and transcriptome analysis reveal recent species-specific gene duplications in the plastic gilthead sea bream (Sparus aurata). Front Mar Sci. 2019;6:760.

  92. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  PubMed  Google Scholar 

  93. Liao Y, Smyth GK, Shi W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 2019;47:e47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinformatics. 2011;27:863–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

    Article  PubMed  PubMed Central  Google Scholar 

  96. Rice P, Longden L, Bleasby A. EMBOSS: The European Molecular Biology Open Software Suite. Trends Genet. 2000;16:276–7.

    Article  CAS  PubMed  Google Scholar 

  97. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  98. Lienhard M, Grimm C, Morkel M, Herwig R, Chavez L. MEDIPS : genome-wide differential coverage analysis of sequencing data derived from DNA enrichment experiments. Bioinformatics. 2014;30:284–6.

    Article  CAS  PubMed  Google Scholar 

  99. Thévenot EA, Roux A, Xu Y, Ezan E, Junot C. Analysis of the human adult urinary metabolome variations with age, body mass index, and gender by implementing a comprehensive workflow for univariate and OPLS statistical analyses. J Prot Res. 2015;14:3322–35.

  100. Weiss S, Van Treuren W, Lozupone C, Faust K, Friedman J, Deng Y, et al. Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME J. 2016;10:1669–81.

    Article  CAS  PubMed Central  Google Scholar 

  101. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  Google Scholar 

  102. Klopfenstein DV, Zhang L, Pedersen BS, Ramírez F, Vesztrocy AW, Naldi A, et al. GOATOOLS: A Python library for Gene Ontology analyses. Sci Rep. 2018;8:10872.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: New features for data integration and network visualization. Bioinformatics. 2011;27:431–2.

    Article  CAS  PubMed  Google Scholar 

  104. Liu W, Xie Y, Ma J, Luo X, Nie P, Zuo Z, et al. IBS: an illustrator for the presentation and visualization of biological sequences. Bioinformatics. 2015;31:3359–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Authors want to thank Dr. Afonso, from Grupo de Investigación en Acuicultura (GIA), IU-ECOAQUA, Universidad de Las Palmas de Gran Canaria, for the PROGENSA breeding program. They acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Funding

Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was supported by the European Union’s Horizon 2020 research and innovation program under grant agreement no. 818367; AquaIMPACT- Genomic and nutritional innovations for genetically superior farmed fish to improve efficiency in European aquaculture. This publication reflects only the authors’ views and the European Union cannot be held responsible for any use that may be made of the information contained herein. Additionally, this study forms part of the ThinkInAzul program and was supported by MCIN with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat Valenciana (THINKINAZUL/2021/024) to JP–S.

Author information

Authors and Affiliations

  1. Nutrigenomics and Fish Growth Endocrinology Group, Institute of Aquaculture Torre de La Sal (IATS, CSIC), 12595, Castellón, Spain

    F. Naya-Català, A. Belenguer, B. Soriano, J. Calduch-Giner & J. Pérez-Sánchez

  2. Grupo de Investigación en Acuicultura (GIA), IU-ECOAQUA, Universidad de Las Palmas de Gran Canaria, Ctra. Taliarte S/N, 35214, Telde, Las Palmas, Canary Islands, Spain

    D. Montero, S. Torrecillas, S. Sarih, M. J. Zamorano & M. Izquierdo

  3. Biotechvana, Parc Científic Universitat de València, 46980, Paterna, Spain

    B. Soriano & C. Llorens

  4. Skretting Aquaculture Research Centre, Stavanger, Norway

    R. Fontanillas

Authors
  1. F. Naya-Català
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. A. Belenguer
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. D. Montero
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. S. Torrecillas
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. B. Soriano
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. J. Calduch-Giner
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. C. Llorens
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. R. Fontanillas
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. S. Sarih
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. M. J. Zamorano
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. M. Izquierdo
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. J. Pérez-Sánchez
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

DM and JP-S planned and oversaw the project. MI contribute with the nutritional programming scheme. RF & DM and MI formulated the diets. ST and SS carried out the experiment. SS, FN-C and JC-G collected and preprocessed tissue samples. The experiment was supervised by DM, MJZ, MI and ST. All analyses were planned and supervised by JP-S. Bioinformatic, genomic and epigenomic analyses were facilitated by FN-C, BS and CL. FN-C, AB and JP-S analyzed and interpreted the data, and generated the figures and tables. FN-C, AB, DM and JP-S wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to J. Pérez-Sánchez.

Ethics declarations

Ethics approval and consent to participate

Animals were bred within the experimental facilities of IU-ECOAQUA (University of Las Palmas de Gran Canaria, Spain) and all procedures of animal rearing and sampling adhered to the ARRIVE guidelines, were approved by the Bioethical Committee of the University of Las Palmas de Gran Canaria (REF: OEBA_ULPGC_26_2019) and were carried out according to European animal directives (2010/63/EU) and Spanish laws (Royal Decree RD53/2013) for the protection of animals used in scientific experiments.

Consent for publication

Not applicable.

Competing interests

RF was employed by the company SKRETTING (Stavanger, Norway). The remaining authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Supplementary Table 1. Detailed sequencing data obtained in this study.

Additional file 2:

Supplementary Table 2. Comparison between RNA-seq results and real-time PCR validation.

Additional file 3:

Supplementary Figure 1. Validation plots of the PLS-DA models, consisting in 500 random permutations, of RNA-seq data for GS (A) and REF (B) animals, and of MBD seq for the 4 experimental groups (C). Heatmap showing the abundance distribution (z-score) of the DE transcripts identified to be driving the separation between dietary groups in GS (D) and REF (E) fish, and heatmap showing the abundance distribution (z-score) of the DM regions identified to be driving the separation between experimental groups (C).

Additional file 4:

Supplementary Table 3A. List of enriched GO-BP obtained from selected DE transcripts (p < 0.05 and VIP ≥ 1) from the PLS-DA using RNAseq data in GS animals.

Additional file 5:

Supplementary Figure 2. Bar chart representing the distribution of DM regions in different genomic areas (intron, exon, promoter and CGI in promoter) within the transcripts showing opposite effects in DNA methylation of their DM regions and their expression in GS (A) and REF (B) fish.

Additional file 6:

Supplementary Table 4. List and annotation of selected DM and discriminant (p < 0.05 and VIP ≥ 1) regions from the 5 supra-categories (classified based on the enriched GO-BP terms) with greater number of genes enriched only in GS fish.

Additional file 7:

Supplementary Figure 3. Organization of the 23 genes selected within the Lipid metabolic process GO-BP term as potential candidate epigenetic markers.

Additional file 8: Supplementary Table 5.

Ingredients and proximate composition of the of the gilthead sea bream broodstock diet.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Naya-Català, F., Belenguer, A., Montero, D. et al. Broodstock nutritional programming differentially affects the hepatic transcriptome and genome-wide DNA methylome of farmed gilthead sea bream (Sparus aurata) depending on genetic background. BMC Genomics 24, 670 (2023). https://doi.org/10.1186/s12864-023-09759-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09759-7

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us