Background

The genetic basis of muscle fat deposition in pigs is not well known. So far, we have only identified a limited number of genes involved in the absorption, transport, storage and catabolism of lipids. Such information is crucial to interpret, from a biological perspective, the results of genome-wide association analyses for intramuscular fat content and composition traits. Herewith, we have investigated how the ingestion of food changes gene expression in the gluteus medius muscle of Duroc pigs.

Results

By comparing the muscle mRNA expression of fasted pigs (T0) with that of pigs sampled 5 h (T1) and 7 h (T2) after food intake, we have detected differential expression (DE) for 148 (T0-T1), 520 (T0-T2) and 135 (T1-T2) genes (q-value <0.05 and a |FC| > of 1.5). Many of these DE genes were transcription factors, suggesting that we have detected the coordinated response of the skeletal muscle to nutrient supply. We also found DE genes with a dual role in oxidative stress and angiogenesis (THBS1, THBS2 and TXNIP), two biological processes that are probably activated in the post-prandial state. Finally, we have identified several loci playing a key role in the modulation of circadian rhythms (ARNTL, PER1, PER2, BHLHE40, NR1D1, SIK1, CIART and CRY2), a result that indicates that the porcine muscle circadian clock is modulated by nutrition.

Conclusion

We have shown that hundreds of genes change their expression in the porcine skeletal muscle in response to nutrient intake. Many of these loci do not have a known metabolic role, a result that suggests that our knowledge about the genetic basis of muscle energy homeostasis is still incomplete.

Physiological genomics aims to understand the molecular basis of highly complex biological processes by applying high-throughput technologies to the large-scale analysis of genomes, transcriptomes and proteomes [1]. We have a very limited understanding of the physiological genomics of intramuscular fat (IMF) content and composition traits in pigs. Several RNA-seq studies comparing the muscle transcriptomes of pigs with divergent lipid profiles have been performed, demonstrating the differential expression of a number of genes related with carbohydrate and lipid metabolism [2–4]. Noteworthy, genome-wide association studies (GWAS) of blood lipid traits in humans have uncovered the existence of a large number of genes strongly associated with plasma lipid concentrations whose involvement in lipoprotein metabolism had never been reported before [5]. For instance, Teslovich et al. [6] performed a GWAS for lipid traits in 100,000 individuals and identified several associated loci (e.g. GALNT2, PPP1R3B, and TTC39B) whose participation in lipid metabolism had not been described previously. Similarly, the Global Lipids Genetics Consortium reported 62 novel loci displaying significant associations with blood lipid levels, and 30 of them had never been previously connected to lipid metabolism [7]. In the light of these results, we can infer that many genes contributing to muscle fat deposition remain to be identified.

The skeletal muscle compartment encompasses a substantial fraction of the body weight and accounts for ≈75% of total insulin-stimulated glucose uptake [8]. Moreover, adipose and muscle tissues absorb most of the chylomicrons generated after a meal consumption [9]. Fat deposition in the porcine muscle may depend, at least in part, on the activation of genes that regulate the uptake, transport, storage, synthesis and degradation of fatty acids (FA) and carbohydrates. As a first step to identify such genes, we have investigated how the profile of pig muscle mRNA expression changes in response to nutrient supply.

Animal material and metabolic profile

A group of 36 female piglets belonging to a commercial Duroc line were brought, after weaning (age = 3–4 weeks), to the IRTA-Pig Experimental Farm at Monells (Girona, Spain). They were fed with a transition feed for 40 days, and, at an approximate age of 2 months, they entered the fattening period. Gilts were housed individually and fed ad libitum with a commercial feeding diet (13% and 5.5% of crude protein and crude fat respectively) until they reached an average live weight of 73 ± 1.2 kg (161 ± 1.1 days). The post-prandial time-points at which muscle gene expression should be analysed were chosen on the basis of the following experiment (experiment 1): we selected at random eight Duroc gilts (out of the 36), with an approximate age of 100 days, and blood samples were taken with citrate Vacutainer tubes before feeding and 2, 4, 6 h. after feeding. These 32 samples were submitted to the Veterinary Clinical Biochemistry Service of the Universitat Autònoma de Barcelona (http://sct.uab.cat/sbcv). The following metabolites were measured using standard protocols: plasma glucose, triglycerides, cholesterol and non-esterified fatty acids.

In experiment 2, we analysed the transcriptomic changes associated with food intake by sequencing the muscle transcriptomes of the 36 Duroc gilts mentioned in the previous paragraph. These gilts were slaughtered at the IRTA-Experimental slaughterhouse in Monells (Girona, Spain) in controlled conditions and complying all national welfare regulations. These 36 sows fasted 12 h prior slaughtering and then 12 of them were stunned, with high concentrations of CO2 to minimize pain, and bled (T0, fasting). The remaining 24 gilts were supplied with a standard feed ad libitum, and slaughtered 5 h (T1, N = 12) and 7 h (T2, N = 12) after T0, following the same procedure reported above. Before slaughter, we took blood samples from these sows and triglyceride and plasma free FA were measured at the Veterinary Clinical Biochemistry Service of the Universitat Autònoma de Barcelona (http://sct.uab.cat/sbcv). After slaughtering, samples of the gluteus medius muscle were collected and submerged in RNAlater (Ambion), being stored at −80 °C until use.

RNA isolation and library construction and sequencing

Each muscle sample was individually submerged in liquid nitrogen and pulverized with a mortar and a pestle. This powder was homogenized with a polytron device in 1 mL of TRI Reagent (Thermo Fisher Scientific, Barcelona, Spain). Total RNA was extracted from gluteus medius muscle samples by using the acid phenol method implemented in the RiboPure kit (Ambion, Austin, TX). Total RNA concentration and purity were assessed with a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Barcelona, Spain), while integrity was checked with a Bioanalyzer-2100 equipment (Agilent Technologies, Inc., Santa Clara, CA). Total RNA samples were submitted to the Centre Nacional d’Anàlisi Genòmica (CNAG, http://www.cnag.cat) for sequencing. Individual libraries for each one of the analysed pigs (N = 36) were prepared using the TruSeq Stranded mRNA Library Preparation Kit (Illumina Inc., CA) according to the protocols recommended by the manufacturer. This level of replication is 4-fold higher than the minimum required (3 individuals/group) in standard RNA-seq studies. Each library was paired-end sequenced (2 × 75 bp) in a HiSeq 2000 platform (Illumina Inc., CA) by using the TruSeq SBS Kit v3-HS (Illumina Inc., CA).

Bioinformatic analyses

Quality control of sequence reads was carried out with the FASTQC software (Babraham Bioinformatics, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We made per-sequence and per-base analyses to filter reads according to the following criteria: sequence-read distribution = 75 bp, 100% coverage in all bases, GC-content ~50%, ~25% of A, T, G and C nucleotide contributions, ambiguous base-content <0.1% and a Phred score higher than 30 (i.e. base-calling accuracy larger than 99.9%). Subsequently, sequences were trimmed for any remaining sequencing adapter by using Trimmomatic v.0.22 [10]. Raw reads were mapped to the pig reference genome (version 10.2-) with the STAR Alignment v.2.5. software [11] by using default parameters and STAR 2-pass alignment steps. The FeatureCounts tool [12] was used to summarize counts of unambiguously mapped reads. The expression of each mRNA was estimated with DESeq2 [13]. This software builds a count matrix Kij (with one row for each gene i and one column for each sample j) encompassing the number of sequencing reads that have been unambiguously mapped to a gene in a sample [13]. The main assumption of this method is that read counts follow a negative binomial distribution with mean μij and dispersion αi [13]. A second important assumption is that genes of similar average expression levels are expected to have a similar dispersion αi value. DeSeq2 calculates final dispersion values by using an empirical Bayes approach that shrinks dispersion estimates towards a set of predicted αi values. When dealing with genes that are poorly expressed, log2 fold-change (FC) estimates can have a high variance due to noisiness issues. To avoid this potential problem, DeSeq2 shrinks log2 fold-change estimates, with an empirical Bayes procedure [13]. Finally, a Wald test is used to infer if shrunken log2 fold-change estimates (and their standard errors) are significantly different from zero. In the Wald test, the shrunken estimate of the log2 fold-change is divided by its standard error, generating a z-statistic that can be compared to a standard normal distribution [13]. Correction for multiple testing is achieved by using a false discovery rate approach [14]. We considered as differentially expressed (DE) those mRNAs displaying a |FC| > 1.5 and a q-value <0.05.

Advaita Bio’s iPathwayGuide (http://www.advaitabio.com/ipathwayguide) and the Cytoscape software [15] combined with the ReactomeFIViz app [16] were used to infer if certain gene ontology terms and pathways are enriched across the sets of DE genes as well as to build biological networks. In order to detect the GO categories that are over- or under-represented in the condition under study, Advaita Bio’s iPathwayGuide uses an impact analysis method that relies on classical statistics but also takes into account other key factors such as the magnitude of each gene’s expression change, their type and position in the given pathways, their interactions, etc. [17]. The ReactomeFIViz application can access the Reactome pathways database in order to do pathway enrichment analysis for a set of genes and visualize hit pathways with the aid of Cytoscape [16]. This application can also access the Reactome Functional Interaction (FI) network to construct a FI sub-network based on a set of genes [16]. In our study, the standard ReactomeFIViz “Gene Set/Mutation Analysis” application was employed to build gene functional interaction networks on the basis of a list of DE genes (q-value <0.05 and a |FC| > of 1.5) and curated pathway information contained in the Reactome database. The functional enrichment analyses for pathways and GO annotations were based on a binomial test [16].

In Experiment 1, measurement of the concentrations of plasma glucose, cholesterol, triglycerides and non-esterified fatty acids revealed that glycaemia and lipidemia peaks took place 2 and 4 h after the 8 Duroc gilts began to eat, a result that was very consistent across individuals (Fig. 1). Eating was also accompanied by a marked decrease of plasma free FA (Fig. 1), a finding that agrees well with the role of these metabolites as a source of energy during fasting. We chose 5 and 7 h post-ingestion as time-points to carry out the analysis of differential expression. Our expectation was that T1 would reflect the process of lipid absorption, while T2 would correspond to a posterior phase in which lipids are stored as triglycerides or catabolized in the β-oxidation pathway to generate ATP. Nevertheless, when we measured the concentrations of triglycerides and plasma free fatty acids in the slaughtered sows forming part of Experiment 2 (Additional file 1: Figure S1), we observed that feeding is associated with an increase in the concentration of triglycerides and a decrease of circulating free FA levels, a result that matches the metabolic profile observed in Experiment 1. However, the kinetics of these two metabolites were not identical to those observed in Experiment 1 because 7 h after feeding triglyceride levels were still peaking. Despite this circumstance, our main comparison (fasting vs fed pigs) remains completely valid.

Fig. 1

Kinetics of the average concentrations of plasma glucose, cholesterol, triglycerides and non-esterified fatty acids (FA) in 8 Duroc pigs at four time points: before eating and 2, 4 and 6 h post-ingestion (p.i)

Results of the Advaita Bio’s iPathwayGuide pathway analysis based on the list of genes that are differentially expressed (q-value <0.05 and |fold-change| > 1.5) in the porcine gluteus medius muscle before (T0) vs 5 h (T1) and 7 h (T2) after eating

T0 vs T1

T0 vs T2

T1 vs T2

Pathway

P-value

Pathway

P-value

Pathway

P-value

Circadian rhythm

1.00E-03

Ribosomea

4.97E-06

Ribosomea

2.84E-13

Circadian entrainment

4.00E-03

Circadian rhythm

8.48E-04

Huntington’s disease

2.84E-04

Cholinergic synapse

4.00E-03

Huntington’s disease

1.00E-03

Parkinson’s disease

7.33E-04

Adrenergic signaling in cardiomyocytes

4.00E-03

Legionellosis

5.00E-03

Oxidative phosphorylationa

8.74E-04

Transcriptional misregulation in cancer

7.00E-03

Parkinson’s disease

6.00E-03

Alzheimer’s disease

1.00E-03

TGF-β signaling pathway

1.30E-02

Viral myocarditis

7.00E-03

Tight junction

1.30E-02

GABAergic synapse

1.50E-02

Malaria

7.00E-03

Metabolic pathwaysa

1.80E-02

Malaria

1.60E-02

p53 signaling pathway

1.00E-02

Herpes simplex infection

1.80E-02

Cardiac muscle contractiona

2.40E-02

Alzheimer’s disease

1.10E-02

p53 signaling pathway

2.50E-02

Herpes simplex infection

2.70E-02

Mineral absorption

1.30E-02

Viral myocarditis

2.90E-02

Fructose and mannose metabolisma

3.20E-02

Toxoplasmosis

1.50E-02

Legionellosis

3.20E-02

Neuroactive ligand-receptor interaction

3.20E-02

PPAR signaling pathway

1.90E-02

Amyotrophic lateral sclerosis (ALS)

3.20E-02

Dopaminergic synapse

3.30E-02

Amyotrophic lateral sclerosis (ALS)

2.20E-02

Sulfur metabolisma

3.60E-02

Alanine, aspartate and glutamate metabolisma

3.50E-02

Sulfur metabolisma

2.40E-02

Arrhythmogenic right ventricular cardiomyopathy (ARVC)

5.00E-02

Glutamatergic synapse

3.60E-02

African trypanosomiasis

2.50E-02

Estrogen signaling pathway

3.70E-02

Transcriptional misregulation in cancer

2.90E-02

Bladder cancer

4.10E-02

Cardiac muscle contractiona

3.30E-02

Dilated cardiomyopathy

4.90E-02

Tight junction

4.90E-02

athe P-value corresponding to the pathway was computed using only over-representation analysis

Post-prandial activation of genes with and without known roles in muscle energy homeostasis

Several of the genes that show the most significant DE between fasted and fed animals (Additional file 2: Table S1) have an established role in metabolism, while for others evidence reported in the literature is more tenuous or even absent. For instance, the 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 (PFKFB3, T0-T1: FC = −3.01, q-value = 1.91E-07) gene can modulate glucose homeostasis by regulating the levels of fructose-2,6-biphosphate [18], and there are substantial evidences that the G0/G1 switch 2 (G0S2, T0-T1: FC = 1.84, q-value = 4.03E-02; T0-T2: FC = 2.06, q-value = 9.35E-04) protein is involved in the regulation of the rate-limiting lipolytic enzyme adipose triglyceride lipase [19].

The analysis of Additional file 2: Table S1 also evidences the existence of DE for several genes with a plausible but poorly characterized role in metabolism. A good example is the mitoguardin 2 (MIGA2, T0-T1: FC = 1.62, q-value = 1.86E-02; T0-T2: FC = 2.22, q-value = 2.10E-05) gene, which shows a dramatic increase in its expression after food intake i.e. MIGA2 is 1.62 and 2.22 times more expressed at 5 and 7 h post-ingestion, respectively. This gene encodes a protein that regulates mitochondrial fusion [20]. Noteworthy, mitochondrial dynamics is highly interconnected with the energy status of the cell, and it has been demonstrated that starvation promotes an acute inhibition of mitochondrial fission [21]. Another gene of interest is syndecan 4 (SDC4, T0-T1: FC = −1.80, q-value = 3.88E-04; T0-T2: FC = −1.82, q-value = 9.59E-04), whose expression levels decreased at 5 h and 7 h after ingestion. In mammals, this gene has been mostly related with cell-matrix adhesion, migration, neuronal development, and inflammation, but studies performed in Drosophila have revealed that it may also have broad effects on the regulation of energy homeostasis [22]. A third example would be the cysteine- serine-rich nuclear protein 1 (CSRNP1, T0-T1: FC = −1.67, q-value = 5.37E-03; T0-T2: FC = −1.75, q-value = 1.07E-02), a molecule that has been mostly related with T-cell immunity [23] and cephalic neural progenitor proliferation [24]. Interestingly, the expression of this molecule is induced by axin, which appears to promote glucose uptake by enhancing the translocation of GLUT4 [25].

The ingestion of food involves changes in the muscle expression of many transcription factors

As shown in Additional file 2: Table S1, we did not detect significant changes in the expression of several genes with a well-established role in lipid uptake (e.g. CD36, lipoprotein lipase), synthesis (e.g. acetyl-CoA carboxylase, fatty acid synthase, diacylglycerol O-acyltransferase 1), transportation (e.g. FA binding proteins) and catabolism (e.g. genes of the β-oxidation pathway). One of the few exceptions to this general trend was the lipase G locus (LIPG, T0-T1: FC = −1.80, q-value = 4.10E-02), which encodes and endothelial lipase modulating lipoprotein metabolism [26]. This gene shows an important drop in its expression levels (1.8 times) 5 h after food intake, a feature that would result in an inhibition of high-density lipoprotein catabolism [26].

Mammalian ribosomes contain 79 different proteins, all of them being encoded by single-copy genes expressed in all tissues [33]. Interestingly, we have detected significant changes in the expression of several ribosomal protein genes (Additional file 2: Table S1). Ribosomal protein genes formed part of the Reactome functional networks shown in Figs. 3 and 4. Moreover, pathways related with ribosomal biogenesis appeared as significant in Table 1 and Additional file 3: Table S2. When nutrients are available, cells tend to activate energy-consuming anabolic pathways whilst under stress or starvation catabolic processes are predominant [33]. Ribosomal biogenesis consumes 60% of cellular energy and this is the key reason why this process is tightly coupled with nutrient supply [34]. The rapamycin (TOR) signalling pathway is deeply involved in coupling ribosome biogenesis with the energy status of the cell by regulating the expression of ribosomal proteins and RNAs [35]. The fundamental role of ribosomal proteins in skeletal muscle metabolism has been illustrated by generating mice where the ribosomal protein S6 cannot be phosphorylated i.e. these mice are viable and fertile but they show muscle weakness and energy deficit [36]. According to our data, these strong changes in the expression of ribosomal protein genes are observed in the T0-T2 and T1-T2 comparisons, but not in T0-T1. Another intriguing observation of our study is that several of these DE ribosomal protein genes are consistently downregulated (e.g. RPS6KA1, RPL35A, RPS23, RPS21, RPL9 and RPL39), a result that is counterintuitive and hard to explain.

Differential expression of genes related with angiogenesis and oxidative stress

The thrombospondin 1 (THBS1, T0-T1: FC = −1.99, q-value = 8.00E-03) and 2 (THBS2, T0-T2: FC = 2.45, q-value = 5.18E-04) and thioredoxin interacting protein (TXNIP, T0-T1: FC = −1.78, q-value = 1.34E-02; T0-T2: FC = −1.79, q-value = 1.13E-02) genes showed significant DE before and after eating (Additional file 2: Table S1). Moreover, they were integrated in the Reactome functional networks depicted in Figs. 2 and 3. These loci have a dual biological role, regulating both angiogenesis and response to oxidative stress. For instance, THBS1 and THBS2 are negative regulators of angiogenesis [37, 38] and their expression is down- and upregulated by oxidative stress, respectively [39, 40]. This feature agrees well with our study, since we found a post-prandial (both at T1 and T2) decreased and increased expression of THBS1 and THBS2, respectively. The TXNIP protein is one of the main regulators of redox homeostasis [41] and also an angiogenic factor [42]. We have observed a diminished expression of this gene after food ingestion, a finding that agrees well with its function as a promoter of oxidative stress and apoptosis [41].

In the mitochondria, oxidative phosphorylation, by which ATP is synthesized as a source of energy, involves the generation of reactive oxygen species (e.g. superoxide, hydrogen peroxide, hydroxyl radical) as a byproduct [43]. This may promote a state of oxidative stress, i.e. an imbalance between oxidants and antioxidants, resulting in cell and tissue damage. Indeed, a single high-fat meal can temporarily impair endothelial function in healthy individuals and this effect is inhibited by antioxidants [44]. Moreover, lipid peroxidation by reactive oxygen species has been suggested as one of the main mechanisms leading to the development of mitochondrial dysfunction and insulin resistance [45]. On the other hand, it is well known that insulin, which is secreted by the pancreas in response to food ingestion, promotes vasodilation and capillary recruitment in the skeletal muscle, an effect mediated by nitric oxide [46]. These actions on the muscle vasculature are fundamental for the maintenance of glucose homeostasis [47]. As a matter of fact, oxidative stress and neovascularization are two tightly linked biological processes i.e. there are evidences that end products of lipid oxidation can bind the Toll-like receptor 2 promoting an angiogenic response [48]. As a whole, DE of THBS1, THBS2 and TXNIP between pre- and post-prandial states probably reflects the combined redox and vascular response of the porcine skeletal muscle to nutrient availability.

A close relationship between nutritional status and the expression of genes integrated in the muscle circadian clock

One of the main results of our experiment was the detection of DE for a set of genes that form part of the peripheral clock that determines the maintenance of circadian rhythms in the skeletal muscle (Figs. 2 and 3, and Additional file 2: Tables S1, Additional file 3: Tables S2 and Additional file 4: Tables S3). Patterns of DE in the two available comparisons (T0-T1 and T0-T2) were consistent i.e. there was an upregulation of ARNTL (T0-T1: FC = 1.87, q-value = 193E-0.4; T0-T2: FC = 2.43, q-value = 2.99E-13) and NR1D1 (T0-T1: FC = 1.61, q-value = 8.30E-03; T0-T2: FC = 1.87, q-value = 9.52E-04), and a downregulation of PER1 (T0-T1: FC = −2.85, q-value = 3.95E-11; T0-T2: FC = −1.83, q-value = 1.12E-0.2), PER2 (T0-T1: FC = −1.67, q-value = 4.33E-04, T0-T2: FC = −2.48, q-value = 7.03E-14), BHLHE40 (T0-T2: FC = −1.77, q-value = 7.87E-0.5), SIK1 (T0-T1: FC = −2.62, q-value = 1.91E-07), CIART (T0-T1: FC = −2.16, q-value = 5.79E-05; T0-T2: FC = −2.35, q-value = 4.52E-06) and CRY2 (T0-T2: FC = −1.60, q-value = 1.28E-0.2). In mammals, the circadian clock is regulated by either the CLOCK-ARNTL or the NPAS2-ARNTL heterodimers depending on the tissue under consideration [49]. These heterodimers activate the transcription of the Period (PER1 and PER2) and Cryptochrome (CRY1 and CRY2) genes [49]. In diurnal species, the PER and CRY complexes accumulate in the cytoplasm during daytime and they are translocated to the nucleus in the evening, thus repressing their own expression through the interaction with CLOCK/ARNTL [49]. The BHLHE40 molecule is a negative regulator of the ARNTL-CLOCK complex [50]. Other clock genes of interest are SIK1, that regulates the entrainment of the circadian clock [51], CIART, whose inactivation increases the circadian period of locomotor activity in mice [52] and NR1D1, a critical regulator of the circadian clock with strong effects on lipid homeostasis [53].

Our data indicate that food ingestion modulates the expression of circadian genes in the porcine skeletal muscle. It might be argued that this DE is just the obvious consequence of slaughtering pigs at different timepoints (T0 = 0 h., T1 = + 5 h. and T2 = + 7 h.). However, studies performed in model species have revealed that the feeding/fasting cycle is one of the main zeitgebers (time cues) synchronizing the skeletal muscle clock [54]. Noteworthy, this clock plays a key role in muscle physiology by regulating the expression of more than one thousand genes mainly involved in metabolic processes [55]. Muscle lipid deposition in pigs could be affected by the expression of these genes because their inactivation in mouse has evidenced numerous metabolic abnormalities including ectopic fat in the muscle, reduced circulating levels of triglycerides and free fatty acids, obesity, hyperlipidemia and severe hepatic steatosis [49]. Besides, SNPs in the human clock genes have been related with abdominal obesity, increase in carbohydrate intake, higher body mass index and metabolic syndrome [56].

Our results indicate that the ingestion of food affects the expression of many transcription factors that are essential for coordinating the metabolic response triggered by the availability of nutrients. Amongst these, clock genes could be particularly important due to their key role in the adequate synchronization of this response as well as because of their broad effects on muscle metabolism. We have also shown that several genes without an evident link with muscle metabolism change their expression in response to nutrient inflow, an observation that suggests that our knowledge about the genetic basis of energy homeostasis in the porcine muscle is still quite limited. Given the close physiological similarity between pigs and humans, data presented in the current study could be also of interest to understand the consequences of food intake on gene expression in this latter species.

Acknowledgements

The authors are indebted to Selección Batallé S.A. for providing the animal material. We gratefully acknowledge to J. Reixach (Selecció Batallé), J. Soler (IRTA), C. Millan (IRTA), A. Quintana (IRTA) and A. Rossell (IRTA) for their collaboration in the experimental protocols and pig management. Thanks also to the CERCA Programme of the Generalitat de Catalunya.

Funding

Part of the research presented in this publication was funded by grants AGL2013-48742-C2-1-R and AGL2013-48742-C2-2-R awarded by the Spanish Ministry of Economy and Competitivity. We also acknowledge the support of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo Ochoa 2016-2019 (SEV-2015-0533) grant awarded to the Center for Research in Agricultural Genomics. We also acknowledge grant 2014 SGR 1528 from the Agency for Management of University and Research Grants of the Generalitat de Catalunya. Tainã F Cardoso was funded with a fellowship from the CAPES Foundation-Coordination of Improvement of Higher Education, Ministry of Education (MEC) of the Federal Government of Brazil. Emilio Mármol-Sánchez was funded with a PhD fellowship FPU15/01733 awarded by the Spanish Ministry of Education and Culture (MECD).

Availability of data and materials

Authors’ contributions

MA and RQ designed the experiment; RQ, JT and MG were responsible for the experimental protocols and generation of animal material; RGP and OG performed the kinetic study; all authors contributed to the obtaining of biological samples; TFC and EM performed RNA extractions; TFC analysed the data; MA and TFC wrote the paper; all authors read and approved the manuscript.

Ethics approval

Animal care, management procedures and blood sampling were performed following national guidelines for the Good Experimental Practices and they were approved by the Ethical Committee of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

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

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.