Exploring Multifunctional Markers of Biological Age in Farmed Gilthead Sea Bream (Sparus aurata). A Transcriptomic and Epigenetic Interplay for an Improved Fish Welfare Assessment Approach

1. Introduction

Aquaculture industry is rapidly expanding and becoming more intensified in order to meet the increased global demand for fish protein [1,2]. This intensification of production may bring about impairments in health and growth performance [3,4,5], which has contributed to increase the ethical concerns of the consumers on aquaculture welfare issues [6,7]. In that sense, it is essential that the monitoring of welfare status might cover a wide-range of welfare indicators, either input-based, referring to the conditions in which the animals are subjected, or outcome-based, which reflect the degree of fulfilment of their welfare needs [8,9,10]. Examples of the first type of indicators are stocking density or water oxygen concentration and temperature, whereas growth performance, survival rates, external appearance (skin/fin erosion), sex ratio or maturation state are outcome-based indicators. Another approach for classifying welfare indicators is to divide them by whether they are practical and easy to use on a farm (operational welfare indicators, OWIs), or complex and difficult without access to a laboratory for evaluation (laboratory welfare indicators, LABWIs). However, the range and diversity of welfare needs are constantly evolving as more knowledge and new technologies become available. In that sense, it is noteworthy that common measurements of cortisol on blood plasma are now moving towards other target tissues (e.g., faeces, scales and mucus) [11] to minimize large variability and dramatic increases by the sampling itself in the aquatic environment [12,13], though each target tissue would provide different information on the welfare state and particularly on the timing of the stress response. Likewise, behavioural monitoring is often difficult in the aquatic environment [14,15], but the use of miniaturized data-loggers externally attached to the operculum of sentinel fish represent a practical solution for a continuous and accurate individual tracking of swimming activity and breathing rates in tank-based rearing systems [16].

The monitoring of the microorganisms living in or around the farmed animals is also a very promising aquaculture welfare indicator that can serve to plan corrective actions to mitigate an environmental issue, or to modulate the physiological state and the response of the organism from a holobiont perspective [17]. Thus, experimental evidence supports the use of some bacteria taxa as strong markers of thermal stress [18,19] or microplastic exposure [20,21] in both salmonid and non-salmonid fish, which could be monitored in a real time and in a cost-effective manner with the advent of the 16S metabarcoding techniques based on the Nanopore technology [22,23], though we are far to establish references values for a wide spectrum of species and culture condition across the production cycle [24,25]. Other welfare indicators that can aid in developing an improved welfare assessment system are those based on age-related cellular and molecular modifications. Certainly, age-associated DNA methylation changes allow the construction of extremely accurate age estimators, termed “epigenetic” or “DNA methylation” clocks, which have been developed for a number of vertebrates, including mammals, birds and fish [26,27,28]. Thus, several piscine DNA methylation clocks have been built, to be applied in fisheries and conservation biology programs, with a focus on searching CpGs that mostly reflect the chronological age [28]. Nevertheless, many age-related epigenetic variations are also sensitive to environmental influences [29], either in mammals [30] or fish [31], and there is accumulating evidence that, depending on the CpG loci used, epigenetic clocks can inform not only about chronological but also biological age [32,33,34], which bring the possibility of a high resolution and precise understanding of the age-related pathology and physiology of an individual [35,36].

According to the above findings, epigenetic markers of biological age have been proposed as reliable markers of cumulative welfare in animal production [37,38,39,40]. However, the integration of data from more than one omics layer will add potential value to the measure [34,41], and the combination of transcriptomic and DNA methylation data may link in a more precise manner age-regulated epigenetic marks with specific cellular functions or fitness traits [41,42]. This also applies to fish, and such integrative omic approach has been used in gilthead sea bream to evaluate the success of nutritional programming in the offspring of broodstock fish fed low fish meal/fish oil diets [43], to reveal genotype-by-temperature interactions in male sex determination and spermatogenesis in zebrafish [44], or to better understand the challenge of moderate hypoxia and temperature increases in Atlantic salmon [45]. Likewise, a targeted transcriptomic/epigenetic approach have disclosed the use of the changing gene expression and DNA methylation pattern of sirtuin 1 as a marker of age- and seasonally-mediated changes of energy metabolism in the skeletal muscle of gilthead sea bream [46]. However, to the best of our knowledge, the integration of wide transcriptomic and DNA methylation data as markers of biological age remains to be largely explored in farmed fish, though experimental evidence highly supported a strong effect of main aquaculture stressors on several biological functions [47,48], through epigenetic marks preferentially located at the gene-promoter regulatory regions [49,50,51]. In agreement with all this, our starting hypothesis was that the integration of wide transcriptomic and epigenomic approaches, filtered by the involvement in multiple functions and differentially methylated CpGs in the gene-promoter regulatory region, would provide a list of putative gene markers of biological age for an extended and precise assessment of cumulative fish welfare. For that purpose, massive gene expression analysis (RNA-seq) combined with a genome-wide DNA methylation (MBD-seq) approach was conducted to disclose the different omics multilayer of white skeletal muscle in one- and three-year old gilthead sea bream, the most highly cultured fish in the Mediterranean region. This allowed to identify up to 10 age-regulated genes, with particular epigenetic marks and multifunction capabilities, as an array of potential biological age gene markers that will require to be further validated in independent studies.

2. Results

2.1. Age-Related Patterns of Gene Expression and Their Discriminant Analysis

Illumina sequencing of mRNA muscle samples from three (S+3) and one (S+1) year old animals generated ~526 million single-end (SE) reads (1 x 75), with an average of ~37.6 million reads per sample (Table S1). After trimming and quality filtering, around 0.4% of all muscle reads were discarded, and the remaining reads ranged between 30.2 million (~2.27 Gb) and 46.9 million (~3.52 Gb) in all samples. Over 90% of these pre-processed reads were mapped against the IATS-CSIC gilthead sea bream reference genome, which retrieved 43,257 expressed coding transcripts (78% of total predicted unique transcripts and corresponding to 20,764 unique descriptions, UD) and 52,844 non-coding transcripts. Differential gene expression analysis between S+3 and S+1 fish resulted in 2,237 transcripts significantly changing between groups when One-way ANOVA (P < 0.05) was applied (Figure 1A), whereas the number of differentially expressed (DE) transcripts decreased to 78 with FDR-adjusted P < 0.05. To better explore the differences due to age, a partial least squares-discriminant analysis (PLS-DA) was conducted with the data filtered by One-Way ANOVA. A clear separation of both groups was observed (Figure 1B), being the segregation mainly due to component 1 (91.6%), in a discriminant model that was validated significantly (explained variance, R2Y=97%, and predicted variance, Q2=88%). Over 1,000 discriminant transcripts (VIP ≥ 1) were identified, being 668 known coding transcripts that were generally down-regulated (75.3 %) in older fish (S+3 vs S+1). Of note, all animals were correctly classified in each group by hierarchical cluster analysis (Figure 1C).

2.2. Age-Related Patterns of Differentially Methylated DNA and Its Discriminant Analysis

Using MBD-seq, ~568 million SE reads (1 x 75), containing at least one methylated CpG, were obtained, with an average of ~40.5 million reads per sample (Table S1). After trimming and quality filtering, around 8% of all muscle reads were discarded, and the remaining reads ranged between 30.2 million (~2.27 Gb) and 46.9 million (~3.52 Gb) in all samples. Over 95% of these pre-processed reads were mapped against the IATS-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.9M 25-bp genomic regions that appeared to contain CG dinucleotides susceptible to be altered by methylation. These methylated regions spanned ~273Mb of the total gilthead sea bream genome, and comprised ~14.1M CpG (77.5% of the total 18.3M CG in the gilthead sea bream genome).

Differential methylation studies detected a total of 107,232 25-bp genomic regions changing (P < 0.05) between S+3 and S+1 groups when One-way ANOVA (P < 0.05) was applied (Figure 2A), while the number of differentially methylated (DM) regions decreased to 164 with FDR-adjusted P < 0.05. The normalized methylation values (rpkm, reads per kilo base per million mapped reads) of these DM regions were used as input in discriminant analysis (PLS-DA) to assess the effect of age over DNA-methylation patterns. Like RNA-seq results, when both groups were compared, the discriminant separation was clear and statistically significant (P < 0.01; Figure 2B). The segregation was mostly due to component 1 (97.9%) and all animals were again correctly classified in each group by the hierarchical clustering analyses (Figure 2C). The resulting model showed percentage values of explained (R2Y) and predicted (Q2) variance around 99% and 97%, respectively (Figure 2B). When VIP threshold (VIP ≥ 1) was applied in the validated PLS-DA model, up to 51,100 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 14,366 DM transcripts (Figure 2B), and most of those regions (58.5 %) presented a hypo-methylated condition with advancing age (S+3 vs S+1).

2.3. Functional Network Analysis of Differentially Expressed Transcripts and of Transcripts with Concomitant Differential Expression and Methylation

When results from S+3 vs S+1 fish were compared and the identified discriminant DE transcripts known to be located in coding areas (668) were matched with those containing DM regions of discriminant value (14,366), only 172 were overlapping (Figure 3A). The over-representation analysis of the list of discriminant DE coding transcripts (668) discerned 442 enriched functions (Gene ontology - Biological process, GO-BP, unique terms; Table S2A). The resulting over-represented functions were clustered in 51 supra-categories (GO-BP ancestors), and the different numbers of DE transcripts within each one are presented in Figure 3B. The GO-BP ancestor Regulation of cellular process showed the largest number of DE transcripts, followed by N compound metabolic process and System development. Interestingly, the expression of most transcripts in all assigned supra-categories (nearly 80 % globally) were down-regulated with age (S+3 vs S+1 animals) and only around 20% presented an up-regulated gene expression.

A functional network analysis of the 172 overlapping transcripts (DE containing DM regions) was also conducted, displaying up to 81 different GO-BP terms (Table S2B). These over-represented functions were clustered in 19 GO-BP ancestors, and the different numbers of DM regions within each one are presented in Figure 3C. Among these regions, around 52 % showed hypo-methylation when S+3 and S+1 fish were compared, and the remaining 48 % were hyper-methylated. This contrasts with the proportions of down- and up-regulated DE transcripts after the enrichment analysis. When the results of the latter analysis (Figure 3B) were compared with those of the enrichment of the overlapping DE transcripts with DM regions (Figure 3C), all GO-BP ancestors identified in the latter appeared in the former, containing most of them a substantial number of DE transcripts.

2.4. The Different Patterns of DNA Methylation and Expression

After the enrichment analysis of the 172 overlapping transcripts, the resulting 108 enriched transcripts were identified and separated in two groups. First, we selected those transcripts with 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). After filtering for transcripts with at least 75 % of their DM regions with a negative correlation, a total of 33 transcripts (33 UD), corresponding to 141 regions, were displayed, and they showed a strong and significant negative relationship between differential methylation and expression (Figure 4A). Regarding in-depth these genes, the position of the DM regions and the number of CpGs are presented in Table S3. Additionally, the genomic organization and number of position of DM CpGs are graphically represented in Figure S1.

The second group contained down-regulated transcripts containing hypo-methylated regions and, again, they were filtered for those with at least 75 % of their DM regions with hypo-methylation. A total of 61 transcripts (57 UD), containing 140 DM regions, were identified. In this case, there was no significant relationship within this group of transcripts between differential methylation and expression (Figure 4B), which indicates that, unlike the first set of genes, no clear pattern on the effect of differential methylation on expression was observed here. The position of the DM regions and the number of CpGs of the latter transcripts are also presented in Table S4, and their genomic organization and number of position of DM CpGs, graphically represented in Figure S2. Overall, most genes (UD) were down-regulated (24 containing hyper-methylated regions and 57 with hypo-methylated CpGs), and only a few presented an up-regulated expression (9 associated to hypo-methylated regions). Only those from the first group, presenting a clearer relationship between gene expression and level of DNA methylation (opposite trend, 24 hyper-methylated and down-regulated and 9 hypo-methylated and up-regulated, Figure 4C), were retained.

2.5. Candidate Markers of Biological Age

The selected genes had been previously assigned, as mentioned, to one of the 19 GO-BP ancestors through the functional network analysis, and the corresponding genes within each biological function are presented in Table 1 (i.e., down-regulated genes linked with hyper-methylated regions or vice versa). The functional processes are ordered by the number of associated genes, being the System development the GO-BP term with the highest number of genes, followed by Anatomical structure development and Anatomical structure morphogenesis. Additionally, these selected genes are ordered by the number of associated functions (Table 2 and Table S3), and the first 14 genes, chosen by their involvement in multiple functions (i.e., at least 4 processes), were selected and are presented in Table 2 (sirt1, NAD-dependent deacetylase sirtuin-1; smad1, Mothers against decapentaplegic homolog 5; col5a1, Collagen alpha-1(V) chain; calcrl, Calcitonin gene-related peptide type 1 receptor-like; nrp2, Neuropilin-2; bmp1, Bone morphogenetic protein 1; ramp1, Receptor activity-modifying protein 1; abi3bp, Target of Nesh-SH3; mical2, Protein MICAL-2; thrb, Thyroid hormone receptor beta; spred2, Sprouty-related, EVH1 domain-containing protein 2; arhgap24, Rho GTPase-activating protein 24; psmd2, 26S proteasome non-ATPase regulatory subunit 2; atp1a2, Sodium/potassium-transporting ATPase subunit alpha-2). Interestingly, among the 14 selected genes, a substantial proportion (71 %; 10 out of 14) contained at least one regulatory DM region (i.e., positioned in promoter, the first intron or the first exon). Thus, those 10 genes with a methylation mark in a regulating region (sirt1, smad1, ramp1, psmd2 – up-regulated; col5a1, calcrl, bmp1, thrb, spred2, atp1a2 – down-regulated) have been proposed as putative candidate markers of biological age.

3. Discussion

Refining approaches to safeguard and optimise the wellbeing of fish can help to increase the profitability and sustainability of aquaculture through the production of more robust fish. There is now a vast array of fish welfare indicators, but there is little consensus on the most adequate criteria to evaluate differences in welfare needs among species, rearing systems or biological features [52,53]. Thus, in the gonochoristic teleost European sea bass, elevated temperatures during early development (thermosensitive period) resulted in a masculinization of the population [54,55]. Conversely, an accelerated male-female sex reversal in the protandrous hermaphrodite gilthead sea bream seems to be the phenotypic expression of a cumulative impaired welfare that is both nutritionally and genetically regulated [56,57]. This latter process exemplifies a biological age deviation as the result of a pseudo-feminization process with a lowered estradiol/11ketosterone plasma ratio in reproductive females, which would denote a failure of the negative feedback inhibition of male-female sex reversal. Likewise, at a cellular and molecular level, the present study disclosed a vast number of age-regulated muscle transcripts by comparing one- and three-year old fish with homogeneous weights within groups and with no sex bias. Indeed, in the hermaphrodite gilthead sea bream nearly all one-year old fish are usually males and in this work the three-year old fish employed for sampling had a similar weight of around 1 kg, as expected for males at that age compared to the heavier females when fed a control diet with high levels of fish meal and fish oil [57]. Overall, this study finally allowed to obtain an enriched list of 10 candidate markers of biological aging after the integration of two omics layers (transcriptomics and epigenetics) and filtering by multifunctional genes with DM CpGs in regulatory regions, which will require a final functional validation in independent studies.

Age-mediated shifts in gene expression have been reported in numerous species [58,59], though the type of genes that are regulated vary among different tissues and organisms [60]. Accordingly, a clear impact of age on the muscle transcriptome was observed in the current study (Figure 1), which is consistent with previous studies in mammals [61,62] and fish [63,64]. Moreover, we observed that most discriminant coding DE transcripts were down-regulated with advancing age, which would reflect the normal declining with age of a wide-range of biological processes [59,65,66], though paradoxically it is not the common feature in a number of previous studies in mammals and fish [61,62,64]. In any case, it is well accepted that aging is related, in mammals [67] and fish (killifish) used a model of vertebrate aging, with an overall hypomethylation that becomes first associated to a down-regulation of DNA methyltransferases. In the present study, this fits well with a global hypo-methylation (58.5 %) among the discriminant DM regions in our age model of one- and three-year old fish, that at the same time would co-exist with a site-specific hypermethylation, as reported upon aging in both humans [69] and other animal models [67,70], including fish [28,71,72].

Methylation of DNA is usually associated with the silencing of gene expression [73]. Nonetheless, as mentioned before, most DE genes presented a lowered expression with age regardless of methylation pattern. Certainly, two different types of transcriptomic and epigenetic interplays were found herein with the age-related down-regulated gene expression. Accordingly, a first set of transcripts presented an opposite trend for DNA methylation and expression (33 transcripts and UD), which means that they showed a down-regulated expression due to hyper-methylation (24 transcripts) or vice versa (9 transcripts), resulting in a significant negative correlation between DNA methylation and their expression level (Figure 4A). Conversely, a second group of 61 down-regulated transcripts (57 UD) were linked to hypo-methylated DM regions. The effect of the differential methylation on gene expression may rely on the position of the epigenetic mark because increased methylated CpG sites in the promoter region of a gene results normally in shutdown of its expression [74,75,76]. Likewise, methylation of the first exon/intron may be also strongly related to the transcriptional silencing [77,78]. Nevertheless, DNA methylation patterns are far more complex than originally thought. For instance, hypo-methylation does not seem to be always an activating epigenetic change, since the loss of DNA methylation in the gene body may be accompanied by formation of repressive chromatin, resulting in a decreased expression [79]. Alternatively, a lower methylation level at gene bodies may result in gene silencing due to nucleosome destabilization in transcribed regions and reduced efficiencies of transcription elongation or splicing [80,81]. This hypo-methylation-induced gene expression shutdown might contribute to a certain degree of instability in the chromatin or nucleosome, which could explain the less clear pattern on the effect of the reduced level of methylation in this second set of retained transcripts (Figure 4B) and suggests that this impact seems to be less reproducible and reliable. Hence, only the transcripts of the first group, whose changes in expression are likely a result of the direct action of enzymes involved in the DNA methylation or demethylation, were kept herein for further selection.

Among the 33 transcripts presenting an opposite trend for DNA methylation and expression, those linked to a greater number of processes (>4 GO-BP ancestors) were retained as a first list of selected markers. Such filtering approach rendered 14 multiple function genes that would be more likely affected by environmental insults (e.g., increased temperatures and hypoxia), which in turn serves to disclose a set of molecular markers that become potentially highly responsive to aging and common aquaculture stressors [47,48,82,83]. In that sense, it must be noted that the resulting list of selected genes are largely involved in developmental processes, which agrees with the fact that aging is viewed as an evolutionary conserved developmental process across mammals [84], and vertebrates in general. Moreover, the observation that most of the above selected genes (10 out of 14; 71%) shared a DM region in their promoter regulatory region further supports a physiological role as age-related markers of environmental perturbations [49,50,51], which makes our final set of 10 genes a reliable list of candidate markers of chronological age deviations.

From a functional point of view, it must be noted that the suppression of cellular senescence by SIRTs is mainly mediated through delaying the age-related telomere attrition, sustaining genome integrity, and promoting DNA damage repair [85,86,87,88]. Certainly, the increasing evidence of SIRTs in the field of aging and age-related diseases indicates that they may provide novel targets for treating diseases associated with aging and extended human lifespan [89]. This becomes especially evident for SIRT1, the most extensively studied mammalian SIRT. In that sense, the observation that the gene sirt1 is included in our final list of biological markers confirms and extend our previous gilthead sea bream study, in which muscle sirt1 was shaped by age at the transcriptional and epigenetic level in a farmed fish trough the production cycle [46]. Moreover, this outcome served as a methodological checkpoint of the procedure that we have used herein for revealing potential markers of chronological age that at the same time have the potential to be regulated at a large extent by a challenging environment. Regarding the gene psmd2 (26S proteasome non-ATPase regulatory subunit 2), it encodes for one of the non-ATPase subunits of the 19S regulator lid within the 26S proteasome, a multicatalytic proteinase complex that is involved in the ubiquitin–proteasome system, recognised as a major intracellular protein degradation system with an important role for muscle homeostasis and health [90]. It has been suggested that proteasome activity decreases with age in several tissues [91,92], including skeletal muscle [93]. Nevertheless, proteasomes might play a relevant role in DNA damage signaling and DNA repair [94], as previously indicated for SIRTs, and this psmd2 gene may be involved in cell cycle and increase cell proliferation [95], which might be related to the high levels of psmd2 significantly associated with age [95] and explain its up-regulation in our model of older fish.

Other two sets of genes with a complementary role with aging are those of bmp1 - smad1 and ramp1- calcrl. First, bmp1 (Bone morphogenetic protein 1) encodes for one of the bone morphogenetic proteins (BMP) that bind to dedicated receptors that in turn phosphorylate BMP-responsive Smad proteins (e.g., Smad1 encoded by gene smad1, mothers against decapentaplegic homolog 5) [96]. BMP signalling is indispensable for the regulation of adult muscle mass in normal and pathological situations [97]. So, it may be speculated that the down-regulation of bmp1 might be related to the progressive reduction in muscle mass with age [98,99], whereas smad1 is up-regulated to maximize the use of the available BMP. The other two linked genes, calcrl (Calcitonin gene-related peptide type 1 receptor-like) and ramp1 (Receptor activity-modifying protein 1), encode for two components of G protein-coupled receptors, specifically the calcitonin receptor-like receptor and a receptor activity-modifying protein, respectively. These receptors mediate the effects of peptide hormones, such as calcitonin gene related peptide (CGRP), which play important roles during adulthood and aging, particularly in the neuromuscular system [100,101]. Indeed, in skeletal muscle, CGRP has been shown to potentiate muscle contraction [102] and increase the rate of blood flow locally following muscle contraction [103,104,105]. In the present study, the expression of calcrl is down-regulated and that of ramp1, up-regulated with age. Thus, it may be proposed that increased transcription of ramp1 might be an effective strategy for increasing the effectiveness of CGRP-induced effects [106], so its expression may be increased with age to compensate the reduction in calcrl transcription. Another down-regulated gene was thrb (Thyroid hormone receptor beta), which encodes for a nuclear hormone receptor for thyroid hormones that participate in skeletal muscle contractile function, myogenesis and muscle regeneration [107]. Thus, its down-regulation may be related to the progressive decline in muscle activity with age [98,99,108].

In our experimental model, another candidate marker that was down-regulated with advancing age is col5a1, (collagen alpha-1(V) chain). Collagen is the main structural protein in the extracellular matrix (ECM) of the skeletal muscle and it is essential for the mechanical support of tissue, although the gene expression of ECM components, including collagenic genes, decreases with age [63,109]. Likewise, the gene spred2 (Sprouty Related EVH1 Domain Containing 2) encodes for one of the proteins that regulate growth factor-induced activation of the mitogen-activated protein kinase (MAPK) cascade [110], which is considered a positive regulator of muscle atrophy [111]. Hence, the lower expression of this gene in older animals, and consequently a lower inhibition of MAPK, could be the cause of an increased muscle atrophy with age. Finally, atp1a2 (Sodium/potassium-transporting ATPase subunit alpha-2) is one of the genes encoding for different Na+, K+ ATPase (NKA) α-isoforms and is mainly expressed in muscle, either in mammals [112,113] or fish [114,115,116]. As a euryhaline teleost, gilthead sea bream may exhibit high NKA activity in seawater [117,118,119,120], and the lower expression of this gene might be linked with the reduction of ATPase activity in muscle with age [121].

In summary, as depicted in Figure 5, an integrative approach of transcriptomics and DNA methylation data from one- and three-year old fish applied over successive selection steps allowed to generate a candidate list of muscle markers of biological age in farmed gilthead sea bream. Previous attempts of age-related markers in fisheries and conservation biology studies have been focused on chronological age markers, but in farmed animals attention is now moved towards indicators of the cumulative burden of environmental experiences, which may accelerate or decelerate the age progression along the production cycle. In comparison to common biomarkers of age, such as telomere length, DNA methylation marks have very often a better predictive ability [122], which would be reinforced herein by the integration of two omics layer (transcriptomics and epigenetics), followed by the use of successive filters to disclose age-related markers highly influenced by cumulative life insults. The result is a reduced list of 10 potential biological age markers that require to be validated in independent studies. Afterwards, the generated knowledge (i.e., the gene expression data from these biomarkers) might be applied through an appropriate estimation algorithm (e.g., multiple linear regression) [123,124] to compute individual’s biological age at a field scale as part of an improved welfare assessment system for auditing farmed fish welfare.

4. Materials and Methods

4.1. Ehics Statement

All procedures of animal rearing and sampling were carried out according to present IATS-CSIC Review Board, European animal directives (2010/63/EU), and Spanish laws (Royal Decree RD53/2013) for the protection of animals used in scientific experiments.

4.2. Experimental Setup

One- (S+1) and three- (S+3) year-old gilthead sea bream (GSB) of Atlantic origin were reared at the indoor experimental facilities of the Institute of Aquaculture Torre de la Sal (IATS-CSIC) in 3000-L tanks under natural photoperiod and temperature conditions at the IATS-CSIC latitude (40°5 N; 0°10E). Water temperature was between 23 and 28 °C in the summer when sampling was conducted, the water oxygen concentration was always higher than 75% saturation, and unionized ammonia remained below toxic levels (< 0.02 mg/L). Fish were fed a standard commercial diet (EFICO YM 568; BioMar, Dueñas, Spain) once a day until visual satiety (5 or 6 times per week depending on fish size). Seven fish per age class (one year old, S+1, 100–115 g body weight; three years old, S+3, 1 kg body weight) were anesthetized with 3-aminobenzoic acid ethyl ester (MS-222, 100 μg/mL), and white skeletal muscle (WSM) was rapidly excised, frozen in liquid nitrogen and stored at − 80 °C until RNA and DNA extraction.

4.3. DNA/RNA Extraction

White skeletal muscle 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 2000cSpectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and DNA integrity was assessed in a 1% agarose gel. Total RNA (70–100 μg) from WSM was extracted with the MagMAXTM-96 Total RNA Isolation Kit (Applied Biosystems, Foster City, CA, United States). The RNA concentration and purity were determined using a Nanodrop 2000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). 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.

4.4. DNA/RNA Illumina Sequencing

For the Methyl-binding domain sequencing (MBD-seq) analysis, 300 ng of DNA from S+1 and S+3 samples were fragmented to 200-550 pb 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 States). Fragmented DNA was then submitted to methylation enrichment using the MethylCollectorTM 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 NextSeq 500 HO sequencer as a 1 × 75 nucleotides single-end (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 PRJNA1096613 (BioSample accession numbers: SAMN40757284-297).

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 NextSeq 500 HO sequencer as 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 PRJNA1096613 (BioSample accession numbers: SAMN40757298-311).

4.5. Bioinformatic 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/). In the case of RNA-seq data, libraries were filtered with Prinseq [125] eliminating those with quality < 18, length < 60 bp, and > 5% of Ns in the sequence. Cleaned reads were mapped against gilthead sea bream reference genome [126] (available at https://seabreamdb.nutrigroup-iats.org/), using TopHat [127]. A representative transcriptome per sample was constructed using Cufflinks, being data quality-checked with CummeRbund [128].

In the case of MBD-seq data, pre-processing was performed with Prinseq [125], eliminating those with quality < 20, 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 [129]. 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 5000 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 [130]. Search parameters for CGI were: length ≥ 200, C+G content ≥ 50%, ratio of observed/expected CpGs ≥ 0.60 and window size = 100.

4.6. Statistical Analysis and Data Filtering

DE transcripts in RNA-seq data were retrieved using DESeq2 at two significance thresholds (P < 0.05 and FDR < 0.05) [131]. DM regions in MBD-seq data were obtained using the MEDIPS R package (P < 0.05) [132] 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 analyses (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 [133]. 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 [134]. Discriminant DE transcripts were matched with those transcripts containing discriminant DM regions (multiomic layer integration) and different bioinformatic analysis and filters (i.e., functional GO enrichment, integrative correlation analysis, multifunction filtering, regulatory region – either promoter, first intron or first exon – filter) were applied over the obtained overlapping transcripts (Figure 6) in order to obtain the list of candidate biological age markers. Over-representation tests of GO-BP terms was implemented in the ShinyGO 0.77 web application [135] and statistical significance was accepted at FDR < 0.05. GO-BP levels and supra-categories were retrieved using GOATools [136]. 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 [137].