Nutrient supply affects the mRNA expression profile of the porcine skeletal muscle

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 CO₂ 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 K _ij (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, log₂ fold-change (FC) estimates can have a high variance due to noisiness issues. To avoid this potential problem, DeSeq2 shrinks log₂ fold-change estimates, with an empirical Bayes procedure [13]. Finally, a Wald test is used to infer if shrunken log₂ fold-change estimates (and their standard errors) are significantly different from zero. In the Wald test, the shrunken estimate of the log₂ 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].

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].

Finally, there is a third category of genes, exemplified by the family with sequence similarity 212, member B (FAM212B, T0-T1: FC = 2.04, q-value = 3.36E-02; T0-T2: FC = 2.68, q-value = 1.13E-06), transmembrane protein 169 (TMEM169, T0-T2: FC = 2.83, q-value = 6.81E-07) and matrix metallopeptidase 25 (MMP25, T0-T2: FC = −2.41, q-value = 7.97E-04) loci, that, to the best of our knowledge, have never been reported to participate in the regulation of energy homeostasis.

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].

We observed DE for many genes encoding transcription factors (Figs. 2 and 3, Additional file 2: Table S1) e.g. the AT-rich interactive domain 5B (ARID5B, T0-T2: FC = −2.31, q-value = 5.98E-04) gene, which influences adipogenesis and also the accumulation of postnatal lipid storage [27]; Kruppel-like factor 5 (KLF5, T0-T2: FC = −1.96, q-value = 1.25E-02), that regulates the expression of genes involved in the β-oxidation of FA [28]; NR4A2, (T0-T1: FC = −2.16, q-value = 8.93E-04), a nuclear orphan receptor that controls the expression of genes related with glucose metabolism [29]; CCAAT/Enhancer Binding Protein δ (CEBPD, T0-T1: FC = −2.33, q-value = 6.37E-05; T0-T2: FC = −1.84, q-value = 1.71E-02) that plays an essential role in adipogenesis [30]; and forkhead box O1 (FOXO1, T0-T1: FC = −1.55, q-value = 2.12E-02; T0-T2: FC = −1.66, q-value = 2.7E-02), which integrates glucose utilization and lipogenesis [31]. In the T0-T2 comparison we found a similar pattern, with DE of genes encoding the nuclear receptor NR4A3 (FC = −2.28, q-value = 1.99E-03), SRY-box 9 (SOX9, FC = −2.28, q-value = 6.84E-05) and BTB and CNC Homology 1, Basic Leucine Zipper (BACH2, FC = −2.45, q-value = 4.61E-05) transcription factors, to mention a few (Figs. 2 and 3, Additional file 2: Table S1). In the T0-T2 comparison (Fig. 3, Additional file 2: Table S1), we also detected an increase in the expression levels of the meteorin (METRNL, FC = 1.77, q-value = 7.33E-03) mRNA that encodes an hormone that promotes energy expenditure and glucose tolerance [32].

Feeding elicits strong changes in the expression of ribosomal protein genes

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].