Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 https://doi.org/10.1186/s40104-019-0412-z RESEARCH Open Access Co-expression network analysis predicts a key role of microRNAs in the adaptation of the porcine skeletal muscle to nutrient supply Emilio Mármol-Sánchez1, Yuliaxis Ramayo-Caldas2, Raquel Quintanilla2, Tainã Figueiredo Cardoso1,3, Rayner González-Prendes4, Joan Tibau2 and Marcel Amills1,5* Abstract Background: The role of non-coding RNAs in the porcine muscle metabolism is poorly understood, with few studies investigating their expression patterns in response to nutrient supply. Therefore, we aimed to investigate the changes in microRNAs (miRNAs), long intergenic non-coding RNAs (lincRNAs) and mRNAs muscle expression before and after food intake. Results: We measured the miRNA, lincRNA and mRNA expression levels in the gluteus medius muscle of 12 gilts in a fasting condition (AL-T0) and 24 gilts fed ad libitum during either 5 h. (AL-T1, N = 12) or 7 h. (AL-T2, N = 12) prior to slaughter. The small RNA fraction was extracted from muscle samples retrieved from the 36 gilts and sequenced, whereas lincRNA and mRNA expression data were already available. In terms of mean and variance, the expression profiles of miRNAs and lincRNAs in the porcine muscle were quite different than those of mRNAs. Food intake induced the differential expression of 149 (AL-T0/AL-T1) and 435 (AL-T0/AL-T2) mRNAs, 6 (AL-T0/AL-T1) and 28 (AL-T0/AL-T2) miRNAs and none lincRNAs, while the number of differentially dispersed genes was much lower. Among the set of differentially expressed miRNAs, we identified ssc-miR-148a-3p, ssc-miR-22-3p and ssc-miR-1, which play key roles in the regulation of glucose and lipid metabolism. Besides, co-expression network analyses revealed several miRNAs that putatively interact with mRNAs playing key metabolic roles and that also showed differential expression before and after feeding. One case example was represented by seven miRNAs (ssc-miR-148a-3p, ssc-miR-151-3p, ssc-miR-30a-3p, ssc-miR-30e-3p, ssc-miR-421-5p, ssc-miR-493-5p and ssc-miR-503) which putatively interact with the PDK4 mRNA, one of the master regulators of glucose utilization and fatty acid oxidation. Conclusions: As a whole, our results evidence that microRNAs are likely to play an important role in the porcine skeletal muscle metabolic adaptation to nutrient availability. Keywords: Co-expression analysis, lincRNAs, microRNAs, Pig, Regulatory impact factor, Skeletal muscle * Correspondence: Marcel.Amills@uab.es 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain 5 Departament de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This 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. Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 2 of 21 Background response to nutrient supply by using mRNA profiles as a The majority of nutrigenomic studies in domestic ani- reference [5]. This analysis took into consideration both mals have investigated the effects of dietary factors on changes in the mean (differential expression) and the the mean expression of messenger RNAs (mRNAs) [1], variance (differential dispersion) of gene expression. whereas the potential consequences of nutrition on the Moreover, we have used a co-expression network ap- expression profiles of microRNAs (miRNAs) and long proach to elucidate potential regulatory interactions be- intergenic non-coding RNAs (lincRNAs) have not been tween expressed miRNAs and differentially expressed explored in depth. Although changes in the expression (DE) mRNA genes as well as to investigate the relation- of porcine genes in response to dietary and genetic fac- ship between gene co-expression modules and meat tors have been reported in previous studies [2–6], the quality and fatty acid composition traits recorded in the regulatory co-expression networks underlying such gluteus medius skeletal muscle of Duroc pigs. changes have not been fully elucidated yet [3, 7, 8]. Moreover, gene expression variance (GEV), also referred Materials and methods as gene dispersion, has been often overlooked, being Animal material and phenotypic recording considered just as experimental noise without any bio- The Duroc pig population used in the current work logical significance [9]. Few methods have been explicitly has been previously described [5]. Thirty-six female designed for modeling GEV across samples in RNA-Seq Duroc piglets were transported to the IRTA-Pig Ex- experiments [10, 11], despite the fact that changes in perimental Farm at Monells (Girona, Spain) after gene expression in response to a specific stimulus might weaning (age = 3–4 weeks). Gilts were kept in transi- have a biologically meaningful individual component tion devices and fed ad libitum with a standard tran- that should not be confounded with experimental noise. sition diet until they reached approximately 2 months Indeed, metabolic responses to nutritional factors are of age (around 18 kg of live weight). Subsequently, all often driven by complex signaling pathways and gene- gilts were transferred to fattening pens, where they to-gene interactions that are not necessarily identical were housed individually and fed ad libitum until across the whole cohort of analyzed biological replicates, reaching approximately 155 d of age. Nutritional de- adding an intrinsic source of variation in gene expres- tails about the feed provided to gilts between 60 and sion patterns that is often ignored or modeled as a con- 155 d have been previously reported in [6]. During stant variable [11]. A widely accepted estimator of GEV fattening (60 to 125 d), gilts received feed ad libitum is the biological coefficient of variation (BCV) [12]. In with 14.6% crude protein, 4.25% crude fat, 4.8% crude contrast with the canonical coefficient of variation (CV), fiber, 4.9% ashes, 0.92% lysine, 0.58% methionine + the BCV effectively integrates both technical and bio- cysteine and 3190 kcal/kg. During the finishing period logical variability, thus avoiding the dependence on (126 to 155 d), gilts were also fed ad libitum with a count size that CV commonly shows. diet containing 14.4% crude protein, 5.53% crude fat, When the expression patterns of two experimental 5.1% crude fiber, 4.9% ashes, 0.86% lysine, 0.53% me- groups are compared, differences in the magnitudes of thionine + cysteine and 3238 kcal/kg. Gilts were average gene expression (differential gene expression) slaughtered in the IRTA Experimental Slaughterhouse and GEV (differential gene dispersion) can be observed. in Monells (Girona, Spain) in accordance with rele- Differential dispersion might be particularly useful to vant Spanish welfare regulations. Before slaughter, the identify regulatory changes induced by the experimental 36 gilts were fasted for 12 h. Subsequently, 12 gilts factor under study. For instance, it is assumed that genes were slaughtered in a fasting condition (AL-T0, N = with low GEV are central members of signal transduc- 12), and the remaining ones were slaughtered 5 h. tion pathways while those with high GEV tend to occupy (AL-T1, N = 12) and 7 h. (AL-T2, N = 12) after re- more peripheral positions in gene networks [13]. How- ceiving food. High concentrations of CO2 were used ever, the central or peripheral position of a given gene in to stun the gilts before bleeding. After slaughter, sam- a network is not necessarily stable across time and it ples of the gluteus medius skeletal muscle were taken could also be altered by the experimental factor being from the 36 gilts, submerged in RNAlater (Thermo analyzed. Differential dispersion could be a useful par- Fisher Scientific, Barcelona, Spain) and stored at − ameter to detect such source of biological variation as 80 °C. The whole experimental design used in the well as to infer its potential consequences. current work is depicted in Fig. 1. In a previous study, we investigated how the patterns Phenotypes listed in Additional file 1: Table S1 were of mRNA expression change in response to food intake recorded in the 36 Duroc gilts. Meat quality traits were by comparing the muscle transcriptomes of fasting vs. measured as described in [14, 15]. Total muscle choles- fed gilts [5]. Herewith, we wanted to determine how the terol content was determined following Cayuela et al. expression profiles of miRNAs and lincRNAs vary in [16], whereas intramuscular fatty acids content and Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 3 of 21 Fig. 1 Depiction of the experimental design used in our study. Gilts were fed ad libitum (N = 36, N = 12 per group) with a commercial feeding diet during the whole growth period. Prior to slaughter, the 36 gilts were fasted for 12 h. The day of slaughter, 12 gilts (AL-T0) were killed under fasting conditions. The remaining 24 gilts were fed during 5 h. (AL-T1) and 7 h. (AL-T2) and they were subsequently slaughtered composition were determined in accordance with previ- RNA > 0.5%) to be sequenced in Sistemas Genómicos ous reports [17]. S.L. (https://www.sistemasgenomicos.com). Individual li- braries for each sample (N = 36) were prepared with the RNA isolation, library preparation and sequencing of TruSeq Small RNA Sample Preparation Kit (Illumina small RNAs Inc., CA) according to the protocols of the manufac- The gluteus medius skeletal muscle RNA-Seq data set turer. Small RNA libraries were then subjected to single- employed in the analysis of lincRNA and mRNA expres- end (1 × 50 bp) sequencing in a HiSeq 2500 platform sion comprised a total of 36 individuals (12 AL-T0, 12 (Illumina Inc., CA). AL-T1 and 12 AL-T2 gilts). Details about the RNA ex- traction and sequencing protocols can be found in [5]. Quality assessment, mapping and count estimation Briefly, gluteus medius skeletal muscle samples were pul- Quality control of paired-end reads was performed with verized and subsequently homogenized in 1 mL of TRI the FASTQC software (Babraham Bioinformatics, http:// Reagent (Thermo Fisher Scientific, Barcelona, Spain). www.bioinformatics.babraham.ac.uk./projects/fastqc/) The RiboPure kit (Ambion, Austin, TX) was used to iso- and filtered reads were trimmed for any remaining se- late the total RNA fraction, and its concentration and quencing adapters with the Trimmomatic v. 0.22 tool purity were determined with a Nanodrop ND-1000 spec- [18], as described in [5, 6]. In the case of single-end se- trophotometer (Thermo Fisher Scientific, Barcelona, quenced reads derived from small RNA molecules, se- Spain). RNA integrity was assessed with a Bioanalyzer- quencing adapters were trimmed and filtered with the 2100 equipment (Agilent Technologies Inc., Santa Clara, Cutadapt software [19], and reads outside a window of CA) by using the Agilent RNA 6000 Nano Kit (Agilent 15–25 nucleotides were discarded. Paired-end trimmed Technologies, Inc., Santa Clara, CA). Libraries were pre- raw reads from RNA-Seq sequences were mapped to the pared with the TruSeq SBS Kit v3-HS (Illumina Inc. CA) porcine Sscrofa.11.1 reference assembly by using the and paired-end sequenced (2 × 75 bp) in a HiSeq 2000 HISAT2 aligner [20] with default parameters. The platform (Illumina Inc., CA) at the Centro Nacional de Stringtie software [21] was subsequently employed to es- Análisis Genómico (https://www.cnag.crg.eu). timate mRNA and lincRNA abundances. Single-end In the present study, we have generated an additional trimmed raw reads derived from small RNAs were also gluteus medius skeletal muscle RNA-Seq data set specif- mapped to the Ssscrofa.11.1 assembly with the Bowtie ically targeting small RNAs and comprising the same 36 Alignment v.1.2.1.1 software [22], and the following individuals cited above. Total RNA was purified as re- specifications for aligning short miRNA reads were taken ported above. The percentage of small-RNA over total into consideration: 1) allowing no mismatches in the RNA was determined with the Agilent Small RNA Kit alignment, 2) removing reads with more than 20 puta- (Agilent Technologies Inc., Santa Clara, CA). All 36 tive mapping sites and 3) reporting first single best samples met the quality threshold (i.e. 0.2–2 μg total stratum alignment (bowtie -n 0 -l 25 -m 20 -k 1 --best RNA with RIN > 7 and miRNA percentage over total --strata). The featureCounts software tool [23] was then Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 4 of 21 used to summarize counts of unambiguously mapped obtained with MDSeq and edgeR were considered to be sig- reads from miRNA-Seq sequences. nificant at a fold change > |1.5| and q-value < 0.05. Differential expression and differential dispersion Gene Ontology and pathway enrichment analysis estimates The lists of mRNA genes detected as DE in the AL-T0/ Raw expression matrices generated on the basis of count AL-T1 and AL-T0/AL-T2 contrasts were used as inputs estimates obtained with Stringtie (mRNAs and lincR- for Gene Ontology (GO) and pathway enrichment ana- NAs) or featureCounts (miRNAs) [21, 23] were normal- lyses. The ClueGO v2.5.0 plug-in application [28] em- ized with the trimmed mean of M-values normalization bedded in the Cytoscape 3.5.1 software [29] was used for method [24]. Sequencing depth and read count per gene determining enriched Reactome and KEGG pathways, as were calculated for each sequenced sample (Add- well as Biological Process enriched GO terms. A two- itional file 15: Figure S1). On the basis of this analysis, sided hypergeometric test of significance was applied for the AL-T0 7197 sample was removed from RNA-Seq determining enriched terms and multiple testing correc- and miRNA-Seq count matrices due to the low read tion for pathway enrichment analyses was implemented coverage observed in the RNA-Seq sequencing data set. with a false discovery rate approach [26], whereas a The presence of influential outliers for each estimate of Bonferroni-based multiple testing correction was used in gene expression was corrected by capping expression the GO enrichment analysis. values laying outside the boundaries of 1.5 times inter- quartile range per gene and fitting them within the 10th Building of co-expression networks and 90th percentiles. For estimating GEV, the BCV was Significant connections between predicted interacting computed for each detected annotated gene as described gene pairs were identified with the Partial Correlation in the edgeR protocol [25], and further discussed in [12]. with Information Theory (PCIT) network inference al- The BCV encapsulates all sources of inter-library vari- gorithm [30]. By using first-order partial correlation ation between replicates, including the contribution of li- coefficients estimated for each trio of genes along brary preparation biases [12]. with an information theory approach, this tool identi- Differentially expressed (DE) and dispersed (DD) genes fies meaningful gene-to-gene putative interactions. were determined by comparing the means and variances of The PCIT approach has been widely used to recon- gene expression in the two AL-T0/AL-T1 and AL-T0/AL- struct co-expression regulatory networks from expres- T2 contrasts. Only mRNAs and miRNAs showing an aver- sion data with good performance [31]. The main aim age expression value above 1 count-per-million (CPM) in at of this analysis was to determine truly informative least 50% (N = 12) of the samples (considering all AL-T0, correlations between node pairs (genes in our con- AL-T1 and AL-T2 samples) were retained for further ana- text), once the influence of other nodes in the net- lyses. Because lincRNAs are much less expressed than work has been considered. mRNAs and miRNAs, all lincRNAs (N = 352) annotated in Pearson pairwise correlation coefficients (r) were the Sscrofa11.1 reference assembly (v. 97) were considered calculated for each expressed miRNA and DE for differential expression and dispersion analyses (a filtering mRNAs in each of the two contrasts (AL-T0/AL-T1 step imposing an expression threshold above 1 CPM would and AL-T0/AL-T2). The pcit function from the have implied the removal of as much as 80% of annotated PCIT R package [30, 32] was then used for detecting lincRNA loci). The edgeR [25] and MDSeq [10] packages meaningful co-expressed gene pairs. To further iden- with default parameters were used for performing differen- tify the putative miRNA-to-mRNA interaction pairs tial expression and dispersion analyses, respectively. The with biological interest, a repressor effect of miRNAs edgeR protocol uses the quantile-adjusted conditional max- on mRNA expression was assumed [33] and, in con- imum likelihood method for detecting differences in gene sequence, only miRNA-to-mRNA co-expressed pairs expression between two groups. Once negative binomial showing r < − 0.5 were retained. Furthermore, we models are fitted to the input counts and dispersion esti- only considered miRNA-to-mRNA interactions with mates are obtained, differential expression is determined by perfect 7mer-m8 pairing between the miRNA-seed using an exact test of significance. Correction for multiple and the 3′-UTR of the putative mRNA targets, hypothesis testing is implemented by using the Benjamini- hence removing spurious miRNA-to-mRNA signifi- Hochberg false discovery rate approach [26]. The MDSeq cant correlations with no robust biological meaning. method implements a re-parametrization of the real-valued To this end, we downloaded the full set of annotated negative binomial distribution to allow the modelling of 3′-UTR sequences in the porcine Sscrofa11.1 assem- gene expression variability [10]. Correction for multiple hy- bly available at BioMart Ensembl repositories (http:// pothesis testing across genes is implemented with the www.ensembl.org/biomart/martview/). Seed portions Benjamini-Yekutieli procedure [27]. The DE and DD genes (2nd to 8th 5′ nucleotides in the mature miRNA) of Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 5 of 21 the annotated set of porcine miRNAs were reverse- values capture different regulatory impact features and complemented and interrogated along the 3′-UTR hence, they can be considered as two independent mea- sequence regions of mRNA genes by making use of surements of the putative relevance of miRNAs as gene the SeqKit toolkit [34]. Additionally, we selected four expression regulators. The RIF1 values for each ith regula- highly expressed and DE miRNAs (ssc-miR-148a-3p, tory factor were calculated as follows: ssc-miR-1, ssc-miR-493-5p and ssc-let-7/ssc-miR-98) and used the TargetScan webserver to evaluate the 1 Xde j¼n evolutionary conservation of their binding sites in RIF1i ¼ PIF j  DW 2ij the 3’UTR of predicted mRNA targets [35]. Only nde j¼1 conserved target mRNAs with TargetScan context++ scores above the 75% percentile were considered as Where nde is the number of DE genes and Phenotype confidently cross-validated. The context++ score de- Impact Factor (PIF) and differential wiring (DW) are de- scribed by Agarwal et al. [35] incorporates the infor- noted by: mation of 14 estimated features in order to rank the 1 2  probability of all the predicted target sites to be bio- PI F j ¼ e1 j −e22j logically functional. 2 DW ij ¼ r1ij −r2ij For those mRNAs predicted to interact with miRNAs, we also investigated if they also interact with other being e1j and e2j the expression of the jth differentially mRNA-encoding genes. In order to focus on relevant expressed gene in both conditions 1 and 2, respectively, putative mRNA-to-mRNA gene interactions, we only whereas r1ij and r2ij represent the co-expression correl- retained those meaningful mRNA co-expressed pairs ation between the ith regulatory factor (miRNAs in our showing |r| > 0.7, as assessed with the PCIT algorithm. case) and the jth DE mRNA gene in conditions 1 and 2, We applied this threshold, which is more stringent than respectively. the one used for miRNA-to-mRNA interactions, because The RIF2 values for each ith regulatory factors were correlations between expressed mRNAs tend to be defined as: higher than those between mRNAs and miRNAs [36]. Hub genes within selected mRNA-to-mRNA gene inter- 1 Xde h 2  2 i j¼n actions (i.e. those mRNAs showing a higher degree of RIF2i ¼ e1 j  r1ij − e2 j  r2ij meaningful connectivity according to the PCIT algo- nde j¼1 rithm), were also identified by calculating a hub score per gene (Ki), defined as: The positive or negative sign of the RIF1 score is mainly determined by the magnitude of the PIF estimates, and hence is dependent on the directionality of the defined xi Ki ¼  contrast (i.e. the AL-T0/AL-T2 vs. AL-T2/AL-T0 con- X trasts would generate RIF1 scores with opposite signs). In Where xi is the number of selected significant connec- contrast, the sign of the RIF2 score reflects the altered tions (|r| > 0.7) reported by the PCIT algorithm and X is ability of the regulators to act as predictors of the abun- the average connectivity within the mRNA-to-mRNA dance of DE genes [37]. co-expression network among DE mRNA genes. Gene co-expression networks were visualized with the Cytos- Association between muscle phenotypes and weighted cape 3.5.1 software [29]. gene co-expression networks Besides, for each selected miRNA-to-mRNA pre- Significant associations between key co-expressed genes dicted interactions, we calculated the Regulatory Im- and meat quality and fatty acids composition traits mea- pact Factor (RIF) of the corresponding miRNAs [37]. sured in the gluteus medius skeletal muscle samples The RIF algorithm aims to identify regulator genes (Additional file 1: Table S1) were determined with the contributing to the observed differential expression in weighted gene correlation network analysis (WGCNA) the analyzed contrasts. Its implementation results in approach [38]. We used the WGCNA R package [38] for two different and inter-connected RIF scores: while building signed weighted gene co-expression modules RIF1 score represents those transcriptional regulators based on mRNA and miRNA genes present in the AL- that are most differentially co-expressed with the T0/AL-T1 and AL-T0/AL-T2 count matrices and dis- most highly abundant and highly DE genes, the RIF2 playing a minimum expression of 1 CPM in at least 50% score highlights those regulators that show the most of samples. Weighted adjacency matrices were built for altered ability to act as predictors of the changes in each expression data set by using a power soft threshold the expression levels of DE genes [37]. Both RIF (β) = 16, as recommended by Langfelder and Horvath Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 6 of 21 [38] for estimating signed correlations based on the assigned to annotated features (genes) after quantification. number of replicates used in our experimental design. With regard to the miRNA-Seq experiment, an average of The obtained weighted adjacency matrices were subse- 8.2 million single-end reads per sample were generated, quently transformed into topological overlapping matri- which were reduced to approximately 6.8 million reads ces (TOM) and corresponding dissimilarities were per sample after quality-check and adapter trimming. calculated to minimize the effect of noise and spurious From these, approximately 77% mapped to the porcine as- co-expression patterns. Hierarchical clustering was then sembly, and an average of 42% single-end mapped reads applied to the dissimilarity matrices (1-TOM) and co- were successfully assigned to annotated microRNAs in the expressed genes were merged into modules through dy- Sscrofa11.1 assembly. The accuracy of the RNA-Seq pro- namic tree branch cutting. Highly inter-connected mod- cedures employed in the current work were previously ules were finally merged by calculating their eigengenes validated by Cardoso et al. [40], analyzing the differential and setting a minimum height cut of 0.25 and a mini- expression of eight genes based on RNA-Seq results and mum module size of 30 genes for each identified gene real-time quantitative PCR measurements of gene expres- co-expression module. sion. Such comparison showed a high concordance be- To further elucidate whether the inferred gene co- tween the results obtained with these two independent expression modules were significantly associated with methods [40]. the variation of meat quality and fatty acids compos- We have characterized and compared the muscle expres- ition traits (Additional file 1: Table S1), module eigen- sion profiles of lincRNAs, miRNAs and mRNAs in three genes (MEs) were defined as the first principal groups of pigs (Fig. 1): AL-T0 (fasted), AL-T1 (5 h after feed- component calculated with the Principal Component ing) and AL-T2 (7 h after feeding). The computed BCVs Analysis (PCA) algorithm. In this way MEs summarize measuring the range of variability in gene expression across the co-expression patterns of all genes within each biological replicates within the same group were markedly module into a single variable. Measured phenotypes elevated for lincRNAs, moderate for mRNAs and low for were then correlated with each defined ME. Correlated miRNAs, which ultimately showed a very stable and homo- phenotype-module pairs were considered to be signifi- geneous expression profile across samples (Fig. 2a). More- cant when P-value < 0.05. Co-expressed miRNA-only over, as expressed by the regularized log2 (Rlog) modules were discarded for further analyses. A Student transformation of gene counts according to Love et al. [41], asymptotic P-value approach was finally used for deter- the average expression of lincRNAs was much lower than mining the significance of the contribution of each gene that of mRNAs, while miRNAs occupied an intermediate within the co-expression modules to the correlation co- position between these two extremes (Fig. 2b). In general, efficient between MEs and each one of the recorded lowly expressed genes displayed higher BCVs than genes phenotypes. Relevant genes within significant modules with high levels of expression (Fig. 3). This pattern was espe- were selected based on the estimates of gene signifi- cially relevant for mRNAs (Fig. 3a), with an average esti- cance (GS, P-value < 0.05) obtained for each phenotype- mated background BCV of 0.53 (i.e. 53% of mean variability module significant association. in gene expression across biological replicates expected for Additionally, hub genes within each detected gene mRNA genes), and lincRNAs (mean BCV = 115%, Fig. 3c). co-expression module showing significant correlations In strong contrast, miRNAs showed a narrow range of gene with phenotypic traits were assessed. WGCNA in- expression variability (mean BCV = 37%). Indeed, we did not ferred networks were converted to edge graphs by detect miRNA genes with extremely high BCV values even using the RNAseqDE wrapper R package (https:// when we considered miRNAs expressed at marginal levels github.com/jtlovell/RNAseqDE). Subsequently, hub below 1 CPM (Fig. 3b). With MDSeq tool [10], we computed scores for each gene in the selected co-expression fold-changes (FC) for dispersion estimates. For each contrast, modules were calculated by computing the scaled log2FC dispersion values were then plotted against log2CPM Kleinberg’s hub centrality score as described in the gene expression values (Fig. 4). In general, protein-coding igraph tool (https://igraph.org) [39]. genes with medium to low expression levels (Fig. 4a) showed higher dispersion FC values than those that were highly Results expressed. This antagonistic relationship was much less obvi- Comparing the expression patterns of coding and non- ous for miRNAs or lincRNAs than for mRNAs (Fig. 4b, c). coding RNAs expressed in the porcine skeletal muscle The RNA-Seq data set employed for mRNA and lincRNA Identification of differentially expressed and dispersed quantification encompassed an average of 48.6 million genes paired-end reads per sample, and approximately 93% of Principal component analysis showed a clear cluster- them mapped successfully to the Sscrofa11.1 assembly. ing of samples according to their group of origin Roughly, 76% of unambiguously mapped reads were (AL-T0, AL-T1 and AL-T2) when we considered Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 7 of 21 Fig. 2 Expression variability and quantification of expression levels of mRNAs, microRNAs and lincRNAs. a Biological Coefficient of Variation (BCV) distribution across transcript types within each analyzed group. b DESeq2 regularized log2 mean expression (rlog) values across transcript types within each analyzed group mRNA expression patterns (Fig. 5a), and this was par- Differential expression and/or dispersion results gener- ticularly true in the AL-T0/AL-T2 contrast. This out- ated with MDSeq and edgeR approaches reflected evi- come agrees well with previously reported results dent changes in the skeletal muscle transcriptomic using the same experimental data [5]. In contrast, the profile of pigs after feed intake. These changes were par- clustering of samples based on their miRNA expres- ticularly intense in the case of mRNA genes, with 149 sion patterns was more diffuse (Fig. 5b), and in the and 435 DE mRNAs in AL-T0/AL-T1 and AL-T0/AL- case of lincRNAs, no evident pattern of clustering T2, respectively (Additional file 2: Table S2). Moreover, was observed (Fig. 5c). This lack of sample clustering 6 and 28 miRNAs (q-value < 0.05; |FC| > 1.5) were classi- could be due, at least in part, to the low and very fied by edgeR as DE in AL-T0/AL-T1 and AL-T0/AL-T2 low numbers of annotated pig miRNAs and lincRNAs, respectively (Table 1), whereas no lincRNAs showed sig- respectively. Moreover, the highly variable expression nificant DE in any of the two contrasts. When we con- of lincRNAs across samples could also contribute to sidered a less stringent FC threshold for miRNAs and this lack of clustering. Joint PCA clustering consider- lincRNAs (|FC| > 1.2), we were able to recover 5 add- ing all three contrast groups is depicted in Add- itional DE miRNAs in the AL-T0/AL-T2 contrast (Table itional file 15: Figure S2. 1). With regard to differential dispersion, 27 and 30 DD As previously said, statistical analyses for DE and DD mRNAs were detected with MDSeq in the AL-T0/AL- miRNA and mRNA genes were restricted to loci with ex- T1 and AL-T0/AL-T2 contrasts, respectively (Add- pression levels above 1 CPM in each contrast and in at itional file 3: Table S3), and several of these mRNAs least 50% (N = 12) of the samples (each contrast includes were also differentially expressed (Additional file 2: 23 samples), whereas all annotated lincRNAs, irrespective Table S2). Few DD miRNAs (i.e. 5 in AL-T0/AL-T1 and of their expression levels, were considered. These filtering 1 in AL-T0/AL-T2) and only two DD lincRNAs (in AL- criteria reduced approximately by half the number of ana- T0/AL-T1) were detected (Table 2). lyzed mRNAs, i.e. 10,648 (AL-T0/AL-T1) and 10,714 (AL- T0/AL-T2,) expressed mRNAs from a total of 22,342 an- Functional annotation and pathway enrichment of notated protein-coding genes were selected for further differentially expressed genes analyses. Regarding miRNAs, 35% of annotated miRNAs A total of 26 Reactome and 8 KEGG significantly did not reach the expression threshold of 1 CPM (286 enriched pathways were detected in the AL-T0/AL-T1 expressed miRNAs out of 442 annotated miRNA genes in contrast, whereas 16 Reactome and 14 KEGG enriched both AL-T0/AL-T1 and AL-T0/AL-T2). pathways were identified for the AL-T0/AL-T2 contrast Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 8 of 21 Fig. 3 Biological Coefficient of Variation (BCV) vs. DESeq2 regularized log2 mean expression (Rlog) of (a) mRNAs, (b) microRNAs and (c) lincRNAs in each of the analyzed groups (AL-T0, AL-T1 and AL-T2) (q-value < 0.05). Gene ontology biological process enrich- signaling. Regarding the GO enriched terms, many bio- ment analyses resulted in 65 and 107 significant GO terms logical processes triggered by nutrient availability after for AL-T0/AL-T1 and AL-T0/AL-T2, respectively. A food intake were activated, such as skeletal muscle differ- complete list of enriched pathways and GO terms is entiation (GO:0035914), carbohydrate biosynthetic shown in Additional files 4: Table S4 (AL-T0/AL-T1) and process (GO:0016051), regulation of gluconeogenesis 5: Table S5 (AL-T0/AL-T2). Among the most highly (GO:0035947), glycogen biosynthetic process (GO: enriched pathways, those related with circadian clock 0005978), gluconeogenesis (GO:0006094), energy reserve regulation appeared in both contrasts, as well as other metabolic process (GO:0006112), activation of transcrip- pathways associated with myogenesis, nuclear receptor tion from RNA polymerase II promoter (GO:0006366), re- transcription or NOTCH1, and interleukin 4 and 13 sponse to lipids (GO:0033993), adipose tissue Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 9 of 21 Fig. 4 Log2 Fold change (FC) of the dispersion values estimated with MDSeq tools vs. log2 mean expression (counts-per-million, CPM) of (a) mRNAs, (b) microRNAs and (c) lincRNAs expression patterns in the AL-T0/AL-T1 (left column) and AL-T0/AL-T2 contrasts (right column) development (GO:006012), regulation of fat cell differenti- of the experimental contrasts. With the PCIT algorithm, ation (GO:0045598), circadian regulation of gene expres- we detected 24 (AL-T0/AL-T1) and 55 (AL-T0/AL-T2) sion (GO:0032922), cellular response to external stimulus miRNAs co-expressed (r < − 0.50) with sets of differentially (GO:0071496), response to starvation (GO:0042594) or expressed putative mRNA targets (Additional file 6: Table regulation of energy homeostasis (GO:2000505), to men- S6). For mRNA-to-mRNA connections, only meaningful tion a few (Additional files 4 and 5: Table S4 and S5). co-expression relationships with |r| > 0.7 were considered (Additional file 7: Table S7). Hub genes showing a high de- Construction of co-expression networks and gree of connectivity were prioritized by means of their esti- measurement of regulatory impact factors mated hub score values (K). A list of selected mRNA genes We also aimed to determine whether the expression of and their K values is available in Additional file 8: Table S8. miRNAs is associated with that of mRNAs in each one Among the genes with the top (5%) hub scores, it is worth Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 10 of 21 Fig. 5 Principal Component Analysis (PCA) clustering of gluteus medius skeletal muscle samples (11 AL-T0, 12 AL-T1 and 12 AL-T2 gilts) according to the expression profiles of (a) mRNAs, (b) microRNAs and (c) lincRNAs mentioning the following ones: (1) AL-T0/AL-T1: Rev-Erb-β hub genes forming part of the co-expression modules (NR1D2), BTB domain and CNC homolog 1 (BACH1), generated with the WGCNA algorithm [38]. By doing ETS proto-oncogene 1 (ETS1) and the cAMP responsive so, we found several genes that in both approaches element binding protein 1 (CREB1), and (2) AL-T0/AL- were identified as top central players in the metabolic T2: secretory carrier membrane protein 2 (SCAMP2), response to food intake. For instance, BACH1 and neuraminidase 3 (NEU3), pyruvate dehydrogenase kinase CREB1 genes were among the top hubs in the 4 (PDK4), fatty acid transport protein 4 (SLC27A4), Blue co-expression module corresponding to the AL- thiamine transporter 1 (SLC19A2), NAD kinase (NADK), T0/AL-T1 contrast (Additional file 9: Table S9). With BTB domain and CNC homolog 2 (BACH2) and ARID respect to AL-T0/ALT2, SCAMP2, NEU3 and PDK4 domain-containing protein 5B (ARID5B). We have also genes within the Green co-expression module were compared the results based on K estimates with the sets of also among the top hub transcripts, whereas BACH2 Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 11 of 21 Table 1 microRNAs detected by edgeR as differentially expressed when comparing AL-T0 (fasted) gilts with their AL-T1 (5 h after eating) and AL-T2 (7 h after eating) counterparts Contrast miRNA log2FCb P-value q-valuec log2CPM AL-T0d log2CPM AL-T1d log2CPM AL-T2d AL-T0/AL-T1a ssc-miR-7-5p 0.9978 7.56E-05 2.16E-02 6.8416 7.5115 – ssc-miR-374a-3p 0.8568 4.73E-04 3.81E-02 6.9201 7.5034 – ssc-miR-7 0.9229 5.26E-04 3.81E-02 6.4206 7.0178 – ssc-miR-148a-3p 0.8989 5.97E-04 3.81E-02 13.8692 14.5105 – ssc-miR-1 0.7686 6.66E-04 3.81E-02 16.8124 17.3183 – ssc-miR-32 1.2420 9.92E-04 4.73E-02 2.9824 3.6098 – AL-T0/AL-T2 a ssc-miR-1285 −2.9830 3.47E-09 9.92E-07 7.9799 – 5.5849 ssc-miR-148a-3p 1.3831 2.39E-06 2.83E-04 13.8692 – 14.9315 ssc-miR-7-5p 1.1592 2.97E-06 2.83E-04 6.8416 – 7.6948 ssc-miR-493-5p 0.7464 3.84E-05 2.37E-03 6.4846 – 7.1191 ssc-miR-7 1.0724 4.14E-05 2.37E-03 6.4206 – 7.1910 ssc-miR-22-3p −0.9814 1.01E-04 4.20E-03 12.6857 – 11.7583 ssc-miR-421-5p 1.2893 1.03E-04 4.20E-03 2.6775 – 3.7359 ssc-miR-758 −0.7536 1.24E-04 4.43E-03 5.4106 – 4.5480 ssc-miR-339 −0.8760 1.68E-04 5.34E-03 2.8919 – 2.0274 ssc-let-7f-1 0.7735 2.36E-04 6.43E-03 14.0981 – 14.7031 ssc-let-7f-5p 0.7641 2.74E-04 6.43E-03 12.4788 – 13.0761 ssc-miR-374a-3p 0.9025 2.75E-04 6.43E-03 6.9201 – 7.5867 ssc-miR-30a-3p 0.6600 3.36E-04 6.43E-03 9.8397 – 10.3833 ssc-miR-151-3p 0.6940 3.37E-04 6.43E-03 12.3832 – 12.9732 ssc-miR-129a-3p −1.3858 4.09E-04 7.05E-03 4.7123 – 3.0830 ssc-miR-296-5p −0.9342 4.79E-04 7.61E-03 5.1094 – 3.9239 ssc-miR-30e-3p 0.6430 7.45E-04 1.12E-02 10.8497 – 11.3840 ssc-miR-98 0.7127 1.24E-03 1.69E-02 9.6818 – 10.2075 ssc-let-7a-1 0.4660 1.13E-03 2.53E-02 13.7761 – 14.1002 ssc-let-7a-2 0.4590 1.43E-03 2.53E-02 12.4875 – 12.8046 ssc-miR-503 0.4912 1.12E-03 2.53E-02 7.8380 – 8.1776 ssc-miR-181c −0.6665 2.02E-03 2.56E-02 3.0770 – 2.3722 ssc-miR-32 1.1189 2.11E-03 2.56E-02 2.9824 – 3.5525 ssc-miR-1 0.6586 2.15E-03 2.56E-02 16.8124 – 17.2479 ssc-miR-450b-3p 0.9689 2.78E-03 2.95E-02 1.3512 – 2.0993 ssc-miR-136-5p 0.9319 2.78E-03 2.95E-02 3.4211 – 3.8968 ssc-miR-7857-3p −1.1003 3.03E-03 3.09E-02 2.3255 – 1.5283 ssc-miR-125b −0.5858 1.88E-03 3.20E-02 13.0432 – 12.3566 ssc-miR-361-5p −0.5109 3.02E-03 4.45E-02 7.4597 – 6.8061 ssc-miR-362 −0.5567 3.45E-03 4.61E-02 6.5327 – 5.8194 ssc-miR-218b 0.7746 4.75E-03 4.62E-02 4.8429 – 5.2796 ssc-miR-532-3p −0.6865 4.84E-03 4.62E-02 7.4422 – 6.5930 ssc-miR-365-3p −0.7367 5.36E-03 4.79E-02 9.9681 – 8.9921 a AL-T0: Duroc gilts in a fasting condition (N = 11); AL-T1: Duroc gilts slaughtered after 5 h of food intake (N = 12); AL-T2: Duroc gilts slaughtered after 7 h of food intake (N = 12) b Log2FC: estimated log2 fold change mean expression levels c q-value: P-value corrected for multiple testing with the Benjamini-Hochberg procedure d Log2CPM: estimated log2 counts-per-million (CPM) mean expression levels in AL-T0, AL-T1 and AL-T2 groups Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 12 of 21 Table 2 microRNAs and lincRNAs detected by MDSeq as differentially expressed when comparing AL-T0 (fasted) gilts with their AL- T1 (5 h after eating) and AL-T2 (7 h after eating) counterparts Contrast log2FCb P-value q-valuec Log2CPM AL-T0d Log2CPM AL-T1d Log2CPM AL-T2d AL-T0/AL-T1a miRNA ssc-miR-17-5p −4.0190 3.20E-06 2.81E-03 6.2654 5.7652 – ssc-miR-186-5p −4.1486 1.66E-06 2.81E-03 8.6343 8.2410 – ssc-miR-362 −3.6875 1.64E-05 9.48E-03 6.5327 5.8885 – ssc-miR-451 −3.6825 2.16E-05 9.48E-03 8.1730 8.1217 – ssc-miR-29a-3p −3.3204 1.16E-04 4.07E-02 9.1335 8.7859 – lincRNA ENSSSCG00000032301 3.3076 3.60E-05 1.32E-02 1.4722 5.8072 – ENSSSCG00000031192 −3.8178 1.59E-04 2.93E-02 5.8864 1.7548 – AL-T0/AL-T2a miRNA ssc-miR-1285 −4.1428 4.60E-06 8.04E-03 7.9799 – 5.5849 a AL-T0: Duroc gilts in a fasting condition (N = 11); AL-T1: Duroc gilts slaughtered after 5 h of food intake (N = 12); AL-T2: Duroc gilts slaughtered after 7 h of food intake (N = 12) b Log2FC: estimated log2 fold change mean dispersion levels c q-value: P-value corrected with the Benjamini-Yekutieli procedure d Log2CPM: estimated log2 counts-per-million (CPM) mean expression levels in AL-T0, AL-T1 and AL-T2 groups and ARID5B occupied intermediate positions in the Besides, after estimating the RIF score for each co- ranking of hub genes (Additional file 9: Table S9). expressed miRNA, results were ranked according to Additionally, we used the TargetScan algorithm to their regulatory relevance. A complete list of all RIF evaluate the accuracy of the miRNA-to-mRNA inter- values for miRNAs is presented in Additional file 12: actions predicted with PCIT and 3′-UTR seed match- Table S12. Moreover, a list of the top 5 ranking posi- ing. Four highly expressed DE miRNAs (ssc-miR- tive and negative regulatory miRNAs according to 148a-3p, ssc-miR-1, ssc-miR-493-5p and ssc-let-7/ssc- their RIF1 and RIF2 scores is presented in Tables 3 miR-98) were selected for this task. From a total of and 4, respectively. Interestingly, we observed a high 30 different mRNA genes predicted to be targets of correspondence between miRNAs classified as DE the selected miRNAs (Additional file 6: Table S6), 14 with the edgeR tool and miRNAs categorized by the showed conserved and putatively valid interactions PCIT and RIF algorithms as meaningful regulators (context++ score > 75% percentile) according to pre- (Tables 1, 3 and 4, Additional files 6 and 12: Table S6 dictions made with the TargetScan algorithm (Add- and S12). For instance, ssc-miR-32, which was DE in itional file 10: Table S10). the two considered contrasts, ranked as the second Particularly interesting was the case of the miRNAs (AL-T0/AL-T1) and third (AL-T0/AL-T2) most rele- predicted to bind the 3′-UTR sequence of the PDK4 vant miRNA in terms of RIF1 (Table 3, mRNA (Additional file 11: Table S11), which happened Additional file 12: Table S12). The DE miRNAs (AL- to be the most highly downregulated gene in the AL-T0/ T0/AL-T2) ssc-miR-339 and ssc-miR-1 were also de- AL-T2 contrast (Additional file 2: Table S2). Among the tected as relevant in terms of RIF1 score (Table 3). 7 predicted miRNAs with putative 7mer-m8 binding When considering RIF2 and AL-T0/AL-T2, the ssc- sites in the PDK4 3′-UTR, only two sites appeared to be miR-1285, ssc-miR-129a-3p, ssc-miR-296-5p, ssc-miR- consistently conserved when compared against the cor- 374a-3p and ssc-miR-7-5p DE miRNAs happened to responding orthologous regions in other phylogenetically be among the top predicted regulators (Table 4). In related species (Additional file 15: Figure S3, Add- the AL-T0/AL-T2 contrast, several additional DE itional file 10: Table S10). Noteworthy, the two con- miRNAs also belonged to the group of the top 10 most served sites are predicted to bind to ssc-miR-148a-3p relevant regulators according to their RIF scores, e.g. and ssc-miR-493-5p, which were two of the most highly ssc-miR-22-3p for RIF1 and ssc-miR-148a-3p or ssc- DE miRNAs in the AL-T0/AL-T2 contrast (Table 1). miR-493-5p for RIF2 (Additional file 12: Table S12). Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 13 of 21 Table 3 Top five positive and negative regulatory microRNAs Table 4 Top five positive and negative regulatory microRNAs according to their Regulatory Impact Factor 1 (RIF1) according to their Regulatory Impact Factor 2 (RIF2) Contrast RIF1 Contrast RIF2 a a AL-T0/AL-T1 AL-T0/AL-T1 ssc-miR-450b-5p 1.7939 ssc-miR-129a-3p 1.8373 ssc-miR-32 1.7041 ssc-miR-219a 1.4996 ssc-miR-136-5p 1.3928 ssc-miR-128 1.4256 ssc-miR-542-3p 1.2969 ssc-miR-503 1.2053 ssc-miR-19a 1.2620 ssc-miR-450b-3p 1.0913 ssc-miR-339-3p −0.9864 ssc-miR-455-5p −0.9408 ssc-miR-421-5p −0.9871 ssc-miR-296-5p −1.0613 ssc-miR-503 −1.1680 ssc-miR-143-3p −1.3923 ssc-miR-326 −1.2569 ssc-miR-542-3p −1.4893 ssc-miR-128 −1.2830 ssc-miR-450b-5p −1.5585 AL-T0/AL-T2a AL-T0/AL-T2a ssc-miR-9858-5p 2.7536 ssc-miR-1285 2.2089 ssc-miR-148b-5p 2.4587 ssc-miR-206 1.7993 ssc-miR-32 2.3825 ssc-let-7d-5p 1.7537 ssc-miR-129a-5p 1.9010 ssc-miR-129a-3p 1.5109 ssc-miR-7139-5p 1.3797 ssc-miR-129a-5p 1.3630 ssc-let-7g −1.0629 ssc-miR-296-5p −1.6368 ssc-miR-130b-5p −1.1300 ssc-miR-374a-3p −1.6758 ssc-miR-339 −1.2069 ssc-miR-148b-5p −1.8280 ssc-miR-1 −1.2630 ssc-miR-7-5p −2.0613 ssc-miR-326 −1.3955 ssc-miR-7139-5p −2.6767 a a AL-T0: Duroc gilts in a fasting condition (N = 11); AL-T1: Duroc gilts AL-T0: Duroc gilts in a fasting condition (N = 11); AL-T1: Duroc gilts slaughtered after 5 h of food intake (N = 12); AL-T2: Duroc gilts slaughtered slaughtered after 5 h of food intake (N = 12); AL-T2: Duroc gilts slaughtered after 7 h of food intake (N = 12) after 7 h of food intake (N = 12) Relationship between weighted gene co-expression modules were: meat redness (a*), pH measured 45 min modules and meat quality and muscle fatty acids post-mortem (PH45GM), linoleic acid content (C18:2- composition traits ω6), arachidonic acid content (C20:4), omega-3 (ω3), The WGCNA algorithm applied to mRNA and miRNA omega-6/omega-3 ratio (ω6/ω3), unsaturated fatty acids expression estimates in the AL-T0/AL-T1 and AL-T0/ content (UFA) and polyunsaturated/saturated fatty acids AL-T2 matrices made possible the identification of 5 ratio (PUFA/SFA) and saturated/unsaturated fatty acids and 10 gene co-expression modules, respectively (Add- ratio (SFA/UFA) (Additional file 14: Table S14). A itional file 15: Figure S4 and S5), excluding miRNA-only detailed list of all analyzed phenotypes is shown in co-expression modules. Among these, the identified Additional file 1: Table S1. P-values measuring the modules for the AL-T0/AL-T1 contrast were signifi- significance of the contribution of each gene within co- cantly associated with the following meat quality and expression modules to significantly correlated pheno- fatty acids composition phenotypes measured in the glu- typic traits can be found in Additional files 13: Table teus medius muscle: meat lightness (L*), intramuscular S13 (AL-T0/AL-T1) and 14: Table S14 (AL-T0/AL-T2). pH (PHGM), intramuscular fat content (GMIMF), pal- mitic acid content (C16:0), linoleic acid content (C18:2- Discussion ω6), arachidonic acid content (C20:4), omega-6 fatty Coding and non-coding RNAs show highly divergent acids content (ω6), omega-6/omega-3 ratio (ω6/ω3), patterns of expression in the porcine muscle polyunsaturated fatty acids content (PUFA) and polyun- By comparing mRNAs, miRNAs and lincRNAs expres- saturated/saturated fatty acids ratio (PUFA/SFA), as sion patterns, we have observed that the expression of shown in Additional file 13: Table S13. Regarding the mRNAs in the porcine skeletal muscle is, on average, AL-T0/AL-T2 contrast, gluteus medius phenotypes substantially higher than that of miRNAs and lincRNAs showing significant associations with co-expression (Fig. 2). This finding was expected because previous Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 14 of 21 studies in humans have reflected the same trend for miRNAs are decoupled to some extent [36]. There are lincRNAs [42, 43] and miRNAs [44]. On the other hand, several factors that explain such circumstance. For in- we have also observed an inverse relationship between stance, miRNAs can be sequestered by pseudogene, the expression means of mRNA and lincRNA genes and mRNA, lincRNA or circular RNA transcripts with re- the magnitude of BCVs (Fig. 3a, c), whereas such trend peated miRNA antisense sequences (the so-called was not obvious for miRNAs (Fig. 3b). miRNA sponges), thus limiting their availability to With regard to differential dispersion, the number of regulate the expression of target RNAs [45–47]. More- DD mRNA and miRNA genes was much lower than that over, compelling evidence has been accumulating dur- of DE mRNA and miRNA genes, indicating that nutrient ing past years highlighting the exceptional stability of supply has a stronger impact on the mean expression of certain miRNAs, which show half-lives of days [48, 49]. genes rather than on their BCV. Of course, these two pa- This long half-life might be explained by the protective rameters are closely related, so decreases in the mean ex- effect of the Argonaute protein in isolating naked pression of genes are usually accompanied by increases in single-stranded small miRNA molecules from exonu- the variance of expression (and vice versa), being such cleases within the cell environment [50]. Besides, miR- trend particularly true for mRNAs and lincRNAs. In con- NAs might localize to cell compartments other than trast, miRNAs showed a very resilient and stable pattern the cytosol, where they exert functions unrelated with of expression across replicates (Figs. 3b and 4b). the modulation of mRNA levels [51]. Last but not least, While nutrient supply induced substantial changes in the expression levels of miRNAs do not necessarily cor- the expression of mRNAs (Additional file 2: Table S2), relate with their functional availability as a part of the the absolute number of DE miRNAs was much lower RNA-induced silencing complex [36]. (Table 1), whereas no DE lincRNAs were detected. This result is probably not due to a limited accuracy of RNA- Differentially expressed and dispersed miRNAs are related Seq in detecting differential gene expression, because with the regulation of key metabolic processes in the previous experiments [40] showed a high consistency be- skeletal muscle tween differential gene expression results obtained with As shown in Tables 1 and 2, several miRNAs were de- RNA-Seq and real time quantitative PCR data in the tected as either being DE and/or DD in the AL-T0/AL- same experimental system. However, it should be taken T1 and AL-T0/AL-T2 contrasts. Among the DE miR- into account that the absolute numbers of annotated NAs, we found that ssc-miR-1 and ssc-miR-148a were porcine miRNAs and lincRNAs are much smaller than two of the most expressed and DE miRNAs in AL-T0/ those of mRNAs. Indeed, when the number of DE genes AL-T1 and AL-T0/AL-T2 contrasts (Table 1), whereas is expressed as a proportion (i.e. number of DE genes/ ssc-miR-7-5p was the most highly differentially upregu- number of total analyzed expressed genes), the total lated miRNA in AL-T1 gilts. Both miR-7 and miR-1 amount of DE mRNAs happened to be 1.39% (AL-T0/ regulate the mTOR-related cell response to nutrient AL-T1) and 4.06% (AL-T0/AL-T2). In the case of miR- availability. For instance, miR-1 was found to be directly NAs, such proportions were 2.09% (AL-T0/AL-T1) and upregulated by the myogenic differentiation 1 (MYOD1) 9.79% (AL-T0/AL-T2). Moreover, the average |FC| of gene [52], which is a transcription factor essential for DE mRNAs was 2.12-fold and 2.02-fold in AL-T0/AL- skeletal muscle development and myocyte fusion [53] T1 and AL-T0/AL-T2 respectively, while for miRNAs, and also functions as a circadian modulator in the per- changes of 1.9-fold (AL-T0/AL-T1) and 1.85-fold (AL- ipheral muscle clock [54]. Noteworthy, MYOD1 was also T0/AL-T2) were detected. In the light of these results, it significantly upregulated in the AL-T0/AL-T2 contrast should be concluded that both mRNAs and miRNAs (Additional file 2: Table S2), a finding that agrees well show consistent patterns of differential expression in re- with the observed upregulation of ssc-miR-1 (Table 1). sponse to food intake, while no conclusive evidence has Additionally, miR-7 has been also associated with the been obtained for lincRNAs. This latter observation Akt-mTOR and PI3K/Akt signaling by targeting the in- could be due to the poor annotation of lincRNAs as well sulin receptor substrate 2 (IRS2) and the phosphoinosi- as to their low expression levels and elevated within tide 3-kinase catalytic subunit δ (PIK3CD) [55, 56], two group expression variability (Figs. 2 and 3c), which ul- genes that are integrated in the coordinated signaling timately would make the differential expression analysis cascade in response to nutrient supply to promote skel- much less powerful to detect significant differences. etal muscle growth and differentiation. Nevertheless, the high variance in the expression of Regarding the miR-148 family, it has been reported lincRNAs contrasted strongly with the stable patterns that these miRNAs play a key role in cholesterol me- of expression across contrasts displayed by miRNAs tabolism [57–59] and insulin homeostasis [60]. In a (Figs. 2 and 3b, c). This high stability might be due to fasting/feeding study resembling ours, Goedeke et al. the fact that the expression and silencing activity of [59] reported that miR-148a binds the 3′-UTR of the Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 15 of 21 low density lipoprotein receptor (LDLR) mRNA lead- miR-148a. The MYF6 gene has been previously associated ing to the accumulation of low-density lipoprotein with the regulation of myogenesis and skeletal muscle cell (LDL) cholesterol in blood plasma. Similar results differentiation [8, 71]. A proliferation modulating function were reported by Rotllan et al. [61]. Furthermore, has also been described for TXNIP [72] as well as for FOSL2 Goedeke et al. [59] suggested that the sterol regula- [73], which is also involved in leptin expression regulation tory element-binding transcription factor 1 (SREBF1) [74], whereas DEPP1 downregulation has been associated may activate the expression of miR-148a by targeting with autophagy inhibition [75]. Moreover, ssc-miR-32 and conserved E-box motifs in the miRNA promoter. In ssc-miR-7-5p, two miRNAs that were differentially upregu- the same study, the role of the ATP-binding cassette lated in AL-T1 gilts (Table 1), were predicted to target sev- 1 (ABCA1) gene in the regulation of high-density eral relevant genes (Additional file 6: Table S6) such as the lipoprotein (HDL) cholesterol levels was explored, and activating transcription factor 3 (ATF3), a key regulator of a binding site for miR-148a in the 3′-UTR of ABCA1 glucose and energy metabolism [76, 77] which was signifi- transcripts was predicted, thus providing a functional cantly downregulated in both AL-T0/AL-T1 and AL-T0/ explanation for the inhibitory effect of miR-148a on AL-T2 contrasts (Additional file 2: Table S2). Other relevant plasma HDL cholesterol levels [59]. Other studies additional transcripts that formed part of the miRNA-to- have also linked miRNAs belonging to the miR-148 mRNA interconnected networks were, to mention a few, the family with angiogenesis and glucose metabolism Kruppel-like factor 15 (KLF15), early growth factor 1 (EGR1) through insulin like growth factor 1 receptor (IGF1R) and ARID domain-containing protein 5B (ARID5B), all of target inhibition [62]. which play key roles in muscle lipid metabolism [8, 78, 79], With respect to other relevant DE miRNAs detected or myogenin (MYOG), a gene that is crucial for muscle de- in our study, the miR-30 family and miR-503 have been velopment and differentiation [80]. described to be involved in skeletal muscle differenti- With regard to AL-T2 gilts, it is worth mentioning the ation and fiber-type composition [63, 64]. Moreover, PDK4 gene, which happened to be the most extremely they also regulate adipogenesis [65], a role that has also downregulated mRNA transcript (Additional file 2: been reported for miR-148a [66] and miR-22 [67]. Fur- Table S2) and was also detected as DD in the AL-T0/ thermore, the observed downregulation of miR-22 after AL-T2 contrast (Additional file 3: Table S3). After food ingestion (Table 1) could be the consequence of reconstructing meaningful miRNA-to-mRNA interac- the active influx of glucose within muscle cells after nu- tions, seven miRNAs (ssc-miR-148a-3p, ssc-miR-151-3p, trient supply. Indeed, the glucose transporter 1 (GLUT1) ssc-miR-30a-3p, ssc-miR-30e-3p, ssc-miR-421-5p, ssc- mRNA is targeted by miR-22 [68]. A similar reasoning miR-493-5p and ssc-miR-503) were predicted to have could be extended to miR-17-5p, which binds to the glu- putative binding sites in the PDK4 3′-UTR (Additional cose transporter 4 (GLUT4) mRNA [69] and that was files 6 and 11: Table S6 and Table S11). Noteworthy, all DD but not DE after feed intake (Table 2). these miRNAs were significantly upregulated in the skel- etal muscle of AL-T2 gilts (Table 1), with the only ex- Relevant miRNA-to-mRNA regulatory interactions in ception of ssc-miR-503, (Table 1). Our findings agree response to nutrient supply well with a cooperative and synergistic interaction be- Co-expression network analyses highlighted that the majority tween the aforementioned miRNAs and the PDK4 of DE miRNAs were also potentially meaningful regulatory mRNA, that would result in its strong downregulation factors (Tables 1, 3 and 4, Additional file 12: Table S12). observed in AL-T2 pigs (Additional file 2: Table S2). Other miRNAs also emerged as potential regulators (Tables 3 Interestingly, among the set of miRNAs significantly co- and 4, Additional file 12: Table S12) despite not being de- expressed with PDK4 mRNAs, and also predicted to tected as significantly DE, a finding that would be in agree- interact with its 3′-UTR, ssc-miR-148a-3p and ssc-miR- ment with the very stable and low expression levels detected 493-5p were two of the most significantly upregulated for most miRNAs (Figs. 2 and 3b). These results evidence miRNAs in AL-T2 gilts (Table 1). Moreover, the Tar- the interest of reconstructing regulatory networks in order to getScan analysis [35] showed that both miRNAs have gain new biological insights that canonical differential expres- evolutionarily conserved binding sites in the 3′-UTR of sion analysis cannot yield [70]. Several critical downregulated the PDK4 gene (Additional file 15: Figure S3, Add- transcription factors in AL-T1 animals were identified as itional file 10: Table S10). We may hypothesize that ssc- potential co-expressed targets of ssc-miR-1 and ssc-miR- miR-148a-3p and ssc-miR-493-5p play a key role in the 148a-3p DE miRNAs (Additional files 2 and 6: Table S2 and downregulation of the PDK4 mRNA after food intake, Table S6), e.g. the myogenic factor 6 (MYF6), FOS-related but such hypothesis still needs experimental verification. antigen 2 (FOSL2) and arrestin domain-containing protein 3 Co-expression network analysis also indicated that the (ARRDC3) for ssc-miR-1, and thioredoxin interacting protein PDK4 gene might interact with a broad array of mRNA (TXNIP) and fasting-induced gene protein (DEPP1) for ssc- transcripts (Fig. 6). Among these, several have been already Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 16 of 21 Fig. 6 Selected miRNA-to-mRNA and mRNA-to-mRNA co-expression network according to the PCIT algorithm in the AL-T0/AL-T2 contrast. Differentially expressed miRNAs and mRNAs were considered. Only significant correlations below − 0.5 for miRNA-to-mRNA and above |0.7| for mRNA-to-mRNA interactions where selected. Red and blue edges indicate negative and positive correlations in the co-expression network, respectively mentioned (MYF6, FOSL2, KLF15, ARID5B, DEPP1, well with the previously described role of BACH2 as a MYOG or TXNIP) while others have not, e.g. aryl hydrocar- transcriptional activator of FOXO1 by binding to its pro- bon receptor nuclear translocator like (ARNTL), forkhead moter region [86–88]. The presence of genes involved in box O1 (FOXO1), neuronal PAS domain protein 2 the maintenance of circadian rhythms (NPAS2, ARNTL (NPAS2), BTB domain and CNC homolog 2 (BACH2) or and PER2) was also relevant, as the expression of the the period circadian regulator 2 (PER2). The PDK4 gene is PDK4 mRNA is subjected to circadian fluctuations in re- one of the master regulators of glucose and lipid metabol- sponse to light shifting and insulin and fatty acids avail- ism [81]. Moreover, the PDK4 protein is located in the ability [89–91]. Noteworthy, the potential implications matrix of the mitochondria and inhibits the pyruvate de- of nutrition in the regulation of the porcine peripheral hydrogenase complex, which catalyzes the conversion of clocks was already discussed in two previous studies pyruvate to acetyl-CoA, and hence it is responsible of the using the very same animal material and experimental decrease in glucose utilization and the upregulation of fatty design reported herewith [5, 40], a result that would be acid oxidation in energy-deprived cells under fasting condi- in agreement with the reconstructed PDK4 miRNA-to- tions [82, 83]. mRNA interaction network reported in this study. The observed coordinated downregulation of both PDK4 and FOXO1 mRNAs in the AL-T0/AL-T2 con- mRNA-to-mRNA hub genes reveal glucose and lipid trast (Additional file 2: Table S2) is consistent with the metabolism changes induced by food intake active energy production and fatty acid synthesis of Hub scoring of meaningful mRNA genes from se- muscle cells in response to nutrient supply, as already lected co-expression interaction networks also allowed reported by Cardoso et al. [5]. In fact, the activation of the identification of several relevant transcripts in- FOXO1 is known to enhance PDK4 transcription by volved in organizing the cell response to nutrient binding to its promoter region [84, 85]. Besides, the availability (Additional file 8: Table S8), and several of BACH2 transcription factor was also predicted to be reg- these were also detected as hub genes in WGCNA ulated by ssc-miR-148a-3p (Additional files 6 and 10: analyses (Additional file 9: Table S9). With respect to Table S6 and Table S10) as well as to interact with both AL-T0/AL-T1, the NR1D2 gene was the most prom- FOXO1 and PDK4 mRNAs (Fig. 6). These findings agree inent hub gene among all other transcripts, despite Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 17 of 21 the fact that it was not detected as DE. This tran- these, the MLXIPL, FASN, SCD, SFRP1, SFRP5 and scription factor and its paralog Rev-Erbα (NR1D1) THRSP genes were also significantly upregulated in contribute to establish links between circadian AL-T2 gilts after feeding (Additional file 2: Table S2). rhythms and cell metabolism regulation [92]. Remark- Interestingly, the active/non-active conformation of ably, other relevant top hub genes were not DE, e.g. the muscle acetylcholine receptor function regulating the BACH1 transcription factor, whose inhibition has motor nerve-muscle communication and muscle con- been associated with an increased protection against traction is tightly associated with the concentration of oxidative stress [93], ETS1, which mediates FOXO1 certain surrounding fatty acid components, contributing acetylation and regulates gluconeogenesis in fasting- to stabilize or destabilize its functionality [101], a feeding cycles [94] or CREB1, an important cofactor phenomenon that could explain the observed association for the peroxisome proliferator-activated receptor γ between its δ subunit (CHRND) and the content of ω-3 coactivator 1-α (PPARGC1A), a gene that plays a key fatty acids and ω6/ω3 content ratio in the gluteus med- role in insulin-mediated glucose uptake [95]. ius, as shown in Additional file 14: Table S14. Regarding hub genes detected in the AL-T0/AL-T2 Other genes that are key regulators of lipid metab- contrast (Additional file 8: Table S8), SCAMP2 has been olism such as SCD, ACACA, FABP4, SFRP1, THRSP related to glucose transporters trafficking during insulin or the hub genes SFRP5 and FASN (Additional file stimulation [96], whereas NEU3, which was also highly 9: Table S9), also clustered in a tight co-expression upregulated in fed gilts (Additional file 2: Table S2), module and they were significantly associated with stimulates insulin sensitivity and glucose tolerance [97]. linoleic and arachidonic fatty acids content in the Other relevant examples are: SLC27A4, responsible for gluteus medius muscle (Additional file 14: Table long chain fatty acids metabolism and trafficking [98], S14). The SFRP5 protein has been thoroughly stud- SLC19A2, also highly downregulated in fed gilts (Add- ied as a central regulator of lipid accumulation and itional file 2: Table S2) and reported as being negatively adipocytes differentiation, which are a result of an regulated by glucose uptake [99], and NADK, a protein increased mitochondrial respiration promoted by that phosphorylates NAD+ to generate NADP+, a metab- SFRP5 blocking of Wnt signaling, hence repressing olite tightly linked with the regulation of circadian Wnt-induced oxidative metabolism [102]. The other rhythms [100]. identified SFRP element (SFRP1) has also been re- These findings agree well with data previously re- ported to be located in a genomic region overlapping ported by Cardoso et al. [5], as well as with enrich- a QTL for meat marbling [103, 104]. Moreover, the ment analyses described in this study (Additional files THRSP, MLXIPL and FASN upregulation detected in 4 and 5: Table S4 and Table S5), where many DE our analyses (Additional file 2: Table S2), as well as genes associated with diverse glucose and lipid metab- their contribution to intramuscular lipid content olism pathways and GO terms were highlighted. (Additional files 9 and 14: Table S9 and Table S14) Other biological processes like muscle proliferation could be a reflection of the intramuscular adipocyte associated to nutrient availability and circadian regu- proliferation triggered by the nutrient supply pro- lation provided compelling evidence about the com- vided to AL-T2 fed gilts [105]. Indeed, the MLXIPL plex machinery triggered in the skeletal muscle to is a key carbohydrate-signaling transcription factor respond to nutrient supply after food ingestion. whose activity is enhanced by glucose metabolites, thus binding to carbohydrate response elements Weighted co-expression analyses revealed hub genes (ChoREs) present in the promoters of several key related with lipids metabolism regulation lipid genes such as FASN [106]. Among the gene co-expression modules detected with the WGCNA approach [38], the so-called Red and Conclusions Purple clusters (Additional file 14: Table S14), In conclusion, we have demonstrated that the profiles corresponding to the AL-T0/AL-T2 contrast, con- of expression of lincRNAs and miRNAs in the gluteus tained several relevant lipid metabolism-related genes medius muscle of pigs are very different than those such as the fatty acid binding protein 4 (FABP4), observed for mRNAs. For instance, the mean and the carbohydrate-responsive element-binding protein variance of gene expression are closely interdependent (MLXIPL), fatty acid synthase (FASN), thyroid hor- parameters in the case of mRNAs, while miRNAs do mone responsive protein (THRSP), stearoyl-CoA desa- not show such trend. We have also demonstrated that turase (SCD), acetyl-CoA carboxylase 1 (ACACA) or feeding induces changes mainly in the mean expres- the secreted frizzled-related proteins 1 and 5 (SFRP1 sion of genes rather than on their expression vari- and SFRP5), as well as other loci such as the cholin- ance, a parameter which remains relatively unaffected ergic receptor nicotinic δ subunit (CHRND). From by nutrient supply. Finally, co-expression network Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 18 of 21 analyses predict that miRNAs and hub mRNA genes expression module association with meat quality and fatty acids composition may play an essential role in the regulation of traits in the AL-T0/AL-T2 contrast as determined with the WGCNA tool. mRNAs showing differential expression upon feeding. Such regulatory interactions predicted with in silico Abbreviations tools should be validated experimentally in order to BCV: Biological Coefficient of Variation; CPM: Counts-Per-Million; CV: Coefficient of Variation; DD: Differentially Dispersed; DE: Differentially verify their occurrence as well as to infer their bio- Expressed; FC: Fold Change; GEV: Gene Expression Variability; GO: Gene logical significance in the context of porcine muscle Ontology; GS: Gene Significance; K: Hub Score; lincRNA: long intergenic non- metabolism and nutrition. coding RNA; MEs: Module Eigengenes; miRNA: microRNA; mRNA: messenger RNA; PCA: Principal Component Analysis; PCIT: Partial Correlation with Information Theory; PIF: Phenotype Impact Factor; r: Pearson Correlation coefficient; RIF: Regulatory Impact Factor; Rlog: Regularized log2; TOM: Topological Overlapping matrix; WGCNA: Weighted Gene Correlation Supplementary information Network Analysis Supplementary information accompanies this paper at https://doi.org/10. 1186/s40104-019-0412-z. Acknowledgements The authors are indebted to Selección Batallé S.A. for providing the animal Additional file 1: Table S1. Meat quality and gluteus medius fatty acids material. We gratefully acknowledge to J. Reixach (Selecció Batallé), J. Soler composition traits recorded in AL-T0, AL-T1 and AL-T2 Duroc gilts. (IRTA), C. Millan (IRTA), A. Quintana (IRTA) and A. Rossell (IRTA) for their collaboration in the experimental protocols and pig management. Thanks Additional file 2: Table S2. Differentially expressed mRNAs in the AL- also to the CERCA Programme of the Generalitat de Catalunya. We would T0/AL-T1 and AL-T0/AL-T2 contrasts. also like to thank dr. Alex Clop for critical comments that greatly improved Additional file 3: Table S3. Differentially dispersed mRNAs in the AL- the content of the manuscript. T0/AL-T1 and AL-T0/AL-T2 contrasts. Additional file 4: Table S4. Pathway Enrichment and GO Enrichment Author’s contributions analyses in the AL-T0/AL-T1 contrast. The authors’ responsibilities were as follows: MA and RQ designed the research. MA, RQ, JT, TFC, RGP and EMS conducted the research. EMS Additional file 5: Table S5. Pathway Enrichment and GO Enrichment analyzed the data. YRC and RGP contributed to the integrative analyses. MA analyses in the AL-T0/AL-T2 contrast. and RQ secured funding for the study. EMS and MA drafted the manuscript. Additional file 6: Table S6. miRNA-to-mRNA significant interactions All authors read and approved the final manuscript. detected with the PCIT algorithm in the AL-T0/AL-T1 and AL-T0/AL-T2 contrasts. Funding Additional file 7: Table S7. mRNA-to-mRNA significant interactions de- The research presented in this publication was funded by grants AGL2013– tected with the PCIT algorithm in the AL-T0/AL-T1 and AL-T0/AL-T2 48742-C2–1-R and AGL2013–48742-C2–2-R awarded by the Spanish Ministry contrasts. of Economy and Competitivity. We also acknowledge the support of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Additional file 8: Table S8. Estimated Hub scores (K) values for genes Severo Ochoa 2016–2019 (SEV-2015-0533) grant awarded to the Centre for differentially expressed in the AL-T0/AL-T1 and AL-T0/AL-T2 contrasts. Research in Agricultural Genomics (CRAG). E. Mármol-Sánchez was funded Additional file 9: Table S9. Scaled Kleinberg’s hub scores per gene for with a PhD fellowship FPU15/01733 awarded by the Spanish Ministry of each WGCNA module significantly correlated with phenotypic traits in Education and Culture (MECD). Y. Ramayo-Caldas is financially supported by the AL-T0/AL-T1 and AL-T0/AL-T2 contrasts. the European Union H2020 Research and Innovation programme under Ma- Additional file 10: Table S10. Identification of conserved mRNA rie Skłodowska-Curie grant (P-Sphere) agreement N° 6655919. T. F. Cardoso targets for selected highly expressed and differentially expressed miRNAs was funded with a fellowship from the CAPES Foundation-Coordination of based on the TargetScan Context++ score. Improvement of Higher Education, Ministry of Education of the Federal Gov- ernment of Brazil. Thanks also to the CERCA Programme of the Generalitat Additional file 11: Table S11. miRNAs predicted to bind the 3′-UTR of de Catalunya. the porcine PDK4 gene. Additional file 12: Table S12. Regulatory Impact Factor scores (RIF1 Availability of data and materials and RIF2) for miRNAs classified by the PCIT algorithm as meaningful The small RNA-Seq data set used and/or analyzed in the current study is avail- regulators in the AL-T0/AL-T1 and AL-T0/AL-T2 contrasts. able at the Sequence Read Archive (SRA) database (BioProject: PRJNA595998). Additional file 13: Table S13. Gene co-expression modules signifi- The previously published RNA-Seq data set was also submitted to the SRA cantly associated with meat quality and fatty acids composition traits ac- database (BioProject: PRJNA386796). cording to the WGCNA algorithm (AL-T0/AL-T1 contrast). Additional file 14: Table S14. Gene co-expression modules signifi- Ethics approval cantly associated with meat quality and fatty acids composition traits ac- Animal care and management procedures followed national guidelines for cording to the WGCNA algorithm (AL-T0/AL-T2 contrast). the Good Experimental Practices and were approved by the Ethical Committee of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA). Additional file 15: Figure S1. Sequencing depth obtained for samples analyzed in each one of the two contrasts (AL-T0/AL-T1 and AL-T0/AL-T2). Fig- Consent for publication ure S2. Joint Principal Component Analysis (PCA) clustering of gluteus medius Not applicable. skeletal muscle samples (11 AL-T0, 12 AL-T1 and 12 AL-T2 samples) according to the expression profiles of (A) mRNAs, (B) microRNAs and (C) lincRNAs. Fig- ure S3. Phylogenetically conserved 7mer-8 m predicted binding sites in the Competing interests 3′- UTR of the pig PDK4 gene for (A) ssc-miR-148a-3p and (B) ssc-miR-493-5p The authors declare that they have no competing interests. porcine miRNAs. The TargetScan software was used for generating conserva- tion graphs across the investigated mammalian species. Red nucleotides show Author details 1 complementary matching base-pairs between the seed of the mature miRNA Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, and the 3’UTR of the pig PDK4 gene. Figure S4. Gene co-expression module Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. 2Animal Breeding and Genetics Program, Institute for Research and Technology in association with meat quality and fatty acids composition traits in the AL-T0/ AL-T1 contrast as determined with the WGCNA tool. Figure S5. Gene co- Food and Agriculture (IRTA), Torre Marimon, 08140 Caldes de Montbui, Spain. 3Present address: Embrapa Pecuária Sudeste, Empresa Brasileira de Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 19 of 21 Pesquisa Agropecuária (EMBRAPA), São Carlos, SP 13560-970, Brazil. 23. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose 4 Department of Animal Science, Universitat de Lleida - Agrotecnio Center, program for assigning sequence reads to genomic features. Bioinformatics. 25198 Lleida, Spain. 5Departament de Ciència Animal i dels Aliments, 2014;30:923–30. Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain. 24. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25. Received: 15 July 2019 Accepted: 4 December 2019 25. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40. 26. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a References practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1. Benítez R, Núñez Y, Óvilo C. Nutrigenomics in farm animals. J Invest 1995;57:289–300. Genomics. 2017;4:1. 27. Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple 2. Puig-Oliveras A, Ramayo-Caldas Y, Corominas J, Estellé J, Pérez-Montarelo D, testing under dependency. Ann Stat. 2001;29:1165–88. Hudson NJ, et al. Differences in muscle transcriptome among pigs 28. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. phenotypically extreme for fatty acid composition. PLoS One. 2014;9:e99720. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene 3. Ayuso M, Fernández A, Núñez Y, Benítez R, Isabel B, Barragán C, et al. ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3. Comparative analysis of muscle transcriptome between pig genotypes 29. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. identifies genes and regulatory mechanisms associated to growth, fatness Cytoscape: a software environment for integrated models of biomolecular and metabolism. PLoS One. 2015;10:e0145162. interaction networks. Genome Res. 2003;13:2498–504. 4. Cardoso TF, Cánovas A, Canela-Xandri O, González-Prendes R, Amills M, 30. Reverter A, Chan EKF. Combining partial correlation and an information Quintanilla R. RNA-seq based detection of differentially expressed genes in the theory approach to the reversed engineering of gene co-expression skeletal muscle of Duroc pigs with distinct lipid profiles. Sci Rep. 2017;7:40005. networks. Bioinformatics. 2008;24:2491–7. 5. Cardoso TF, Quintanilla R, Tibau J, Gil M, Mármol-Sánchez E, González- 31. Bellot P, Olsen C, Salembier P, Oliveras-Vergés A, Meyer PE. NetBenchmark: a Rodríguez O, et al. Nutrient supply affects the mRNA expression profile of bioconductor package for reproducible benchmarks of gene regulatory the porcine skeletal muscle. BMC Genomics. 2017;18:603. network inference. BMC Bioinformatics. 2015;16:312. 6. Ballester M, Amills M, González-Rodríguez O, Cardoso TF, Pascual M, 32. Watson-Haigh NS, Kadarmideen HN, Reverter A. PCIT: an R package for González-Prendes R, et al. Role of AMPK signalling pathway during weighted gene co-expression networks based on partial correlation and compensatory growth in pigs. BMC Genomics. 2018;19:682. information theory approaches. Bioinformatics. 2010;26:411–3. 7. Jia C, Kong X, Koltes JE, Gou X, Yang S, Yan D, et al. Gene co-expression network analysis unraveling transcriptional regulation of high-altitude 33. Bartel DP. Metazoan MicroRNAs. Cell. 2018;173:20–51. adaptation of Tibetan pig. PLoS One. 2016;11:e0168161. 34. Shen W, Le S, Li Y, Hu F. SeqKit: a cross-platform and ultrafast toolkit for 8. Muñoz M, García-Casco JM, Caraballo C, Fernández-Barroso MA, Sánchez- FASTA/Q file manipulation. PLoS One. 2016;11:e0163962. Esquiliche F, Gómez F, et al. Identification of candidate genes and 35. Agarwal V, Bell GW, Nam J-W, Bartel DP. Predicting effective microRNA regulatory factors underlying intramuscular fat content through longissimus target sites in mammalian mRNAs. eLife. 2015;4:e05005. dorsi transcriptome analyses in heavy Iberian pigs. Front Genet. 2018;9:608. 36. Mayya VK, Duchaine TF. On the availability of microRNA-induced silencing 9. Komurov K, Ram PT. Patterns of human gene expression variance show strong complexes, saturation of microRNA-binding sites and stoichiometry. Nucleic associations with signaling network hierarchy. BMC Syst Biol. 2010;4:154. Acids Res. 2015;43:7556–65. 10. Ran D, Daye ZJ. Gene expression variability and the analysis of large-scale 37. Reverter A, Hudson NJ, Nagaraj SH, Pérez-Enciso M, Dalrymple BP. RNA-seq studies with the MDSeq. Nucleic Acids Res. 2017;45:e127. Regulatory impact factors: unraveling the transcriptional regulation of 11. Ma C, Ji T. Detecting differentially expressed genes for syndromes by complex traits from expression data. Bioinformatics. 2010;26:896–904. considering change in mean and dispersion simultaneously. BMC 38. Langfelder P, Horvath S. WGCNA: and R package for weighted correlation Bioinformatics. 2018;19:330. network analysis. BMC Bioinformatics. 2008;9:559. 12. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of 39. Csardi G, Nepusz T. The igraph software package for complex network multifactor RNA-Seq experiments with respect to biological variation. research. InterJournal Complex Syst. 2006;1695. http://igraph.org. Nucleic Acids Res. 2012;40:4288–97. 40. Cardoso TF, Quintanilla R, Castelló A, Mármol-Sánchez E, Ballester M, 13. Chalancon G, Ravarani C, Balaji S, Martinez-Arias A, Aravind L, Jothi R, et al. Jordana J, et al. Analysing the expression of eight clock genes in five tissues Interplay between gene expression noise and regulatory network from fasting and fed sows. Front Genet. 2018;9:475. architecture. Trends Genet. 2012;28:221–32. 41. Love MI, Huber W, Anders S. Moderated estimation of fold change and 14. Eusebi PG, González-Prendes R, Quintanilla R, Tibau J, Cardoso TF, Clop A, dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. et al. A genome-wide association analysis for carcass traits in a commercial 42. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, et al. Duroc pig population. Anim Genet. 2017;48:466–9. Integrative annotation of human large intergenic noncoding RNAs reveals 15. Mármol-Sánchez E, Quintanilla R, Jordana J, Amills M. An association analysis for 14 global properties and specific subclasses. Genes Dev. 2011;25:1915–27. candidate genes mapping to meat quality quantitative trait loci in a Duroc pig 43. Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, et al. The landscape of population reveals that the ATP1A2 genotype is highly associated with muscle long noncoding RNAs in the human transcriptome. Nat Genet. 2015;47:199–208. electric conductivity. Anim Genet. 2019. https://doi.org/10.1111/age.12864. 44. de Rie D, Abugessaisa I, Alam T, Arner E, Arner P, Ashoor H, et al. An 16. Cayuela JM, Garrido MD, Bañón SJ, Ros JM. Simultaneous HPLC analysis of α- integrated expression atlas of miRNAs and their promoters in human and tocopherol and cholesterol in fresh pig meat. J Agric Food Chem. 2003;51:1120–4. mouse. Nat Biotechnol. 2017;35:872–8. 17. Mach N, Devant M, Díaz I, Font-Furnols M, Oliver MA, García JA, et al. 45. Ebert MS, Sharp PA. Emerging roles for natural microRNA sponges. Curr Biol. Increasing the amount of n-3 fatty acid in meat from young Holstein bulls 2010;20:R858–61. through nutrition. J Anim Sci. 2006;84:3039–48. 46. Migault M, Donnou-Fournet E, Galibert M-D, Gilot D. Definition and 18. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina identification of small RNA sponges: focus on miRNA sequestration. sequence data. Bioinformatics. 2014;30:2114–20. Methods. 2017;117:35–47. 19. Martin M. Cutadapt removes adapter sequences from high-throughput 47. Pan X, Wenzel A, Jensen LJ, Gorodkin J. Genome-wide identification of sequencing reads. EMBnet.journal. 2011;17:10. clusters of predicted microRNA binding sites as microRNA sponge 20. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low candidates. PLoS One. 2018;13:e0202369. memory requirements. Nat Methods. 2015;12:357–60. 48. Bail S, Swerdel M, Liu H, Jiao X, Goff LA, Hart RP, et al. Differential regulation 21. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. of microRNA stability. RNA. 2010;16:1032–9. StringTie enables improved reconstruction of a transcriptome from RNA-seq 49. Guo Y, Liu J, Elfenbein SJ, Ma Y, Zhong M, Qiu C, et al. Characterization of the reads. Nat Biotechnol. 2015;33:290–5. mammalian miRNA turnover landscape. Nucleic Acids Res. 2015;43:2326–41. 22. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient 50. Lv Z, Wei Y, Wang D, Zhang C-Y, Zen K, Li L. Argonaute 2 in cell-secreted alignment of short DNA sequences to the human genome. Genome Biol. microvesicles guides the function of secreted miRNAs in recipient cells. 2009;10:R25. PLoS One. 2014;9:e103599. Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 20 of 21 51. Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat Rev 75. Salcher S, Hermann M, Kiechl-Kohlendorfer U, Ausserlechner MJ, Obexer P. Mol Cell Biol. 2019;20:21–37. C10ORF10/DEPP-mediated ROS accumulation is a critical modulator of 52. Sun Y, Ge Y, Drnevich J, Zhao Y, Band M, Chen J. Mammalian target of FOXO3-induced autophagy. Mol Cancer. 2017;16:95. rapamycin regulates miRNA-1 and follistatin in skeletal myogenesis. J Cell 76. Lee Y-S, Sasaki T, Kobayashi M, Kikuchi O, Kim H-J, Yokota-Hashimoto H, Biol. 2010;189:1157–69. et al. Hypothalamic ATF3 is involved in regulating glucose and energy 53. Zhang Y, Yu B, He J, Chen D. From nutrient to microRNA: a novel insight metabolism in mice. Diabetologia. 2013;56:1383–93. into cell signaling involved in skeletal muscle development and disease. Int 77. Allison MB, Pan W, MacKenzie A, Patterson C, Shah K, Barnes T, et al. J Biol Sci. 2016;12:1247. Defining the transcriptional targets of leptin reveals a role for Atf3 in leptin 54. Hodge BA, Zhang X, Gutierrez-Monreal MA, Cao Y, Hammers DW, Yao Z, action. Diabetes. 2018;67:1093–104. et al. MYOD1 functions as a clock amplifier as well as a critical co-factor for 78. Boyle KB, Hadaschik D, Virtue S, Cawthorn WP, Ridley SH, O’Rahilly S, et al. downstream circadian gene expression in muscle. eLife. 2019;8:e43017. The transcription factors Egr1 and Egr2 have opposing influences on 55. Kefas B, Godlewski J, Comeau L, Li Y, Abounader R, Hawkinson M, et al. adipocyte differentiation. Cell Death Differ. 2009;16:782–9. microRNA-7 inhibits the epidermal growth factor receptor and the Akt 79. Haldar SM, Jeyaraj D, Anand P, Zhu H, Lu Y, Prosdocimo DA, et al. Kruppel- pathway and is down-regulated in glioblastoma. Cancer Res. 2008;68:3566–72. like factor 15 regulates skeletal muscle lipid flux and exercise adaptation. 56. Fang Y, Xue J-L, Shen Q, Chen J, Tian L. MicroRNA-7 inhibits tumor growth Proc Natl Acad Sci U S A. 2012;109:6739–44. and metastasis by targeting the phosphoinositide 3-kinase/Akt pathway in 80. Ganassi M, Badodi S, Ortuste Quiroga HP, Zammit PS, Hinits Y, Hughes SM. hepatocellular carcinoma. Hepatology. 2012;55:1852–62. Myogenin promotes myocyte fusion to balance fibre number and size. Nat 57. Goedeke L, Aranda JF, Fernández-Hernando C. microRNA regulation of Commun. 2018;9:4232. lipoprotein metabolism. Curr Opin Lipidol. 2014;25:282–8. 81. Jeong JY, Jeoung NH, Park K-G, Lee I-K. Transcriptional regulation of 58. Wagschal A, Najafi-Shoushtari SH, Wang L, Goedeke L, Sinha S, deLemos AS, pyruvate dehydrogenase kinase. Diabetes Metab J. 2012;36:328–35. et al. Genome-wide identification of microRNAs regulating cholesterol and 82. Holness MJ, Sugden MC. Regulation of pyruvate dehydrogenase complex triglyceride homeostasis. Nat Med. 2015;21:1290–7. activity by reversible phosphorylation. Biochem Soc Trans. 2003;31:1143–51. 59. Goedeke L, Rotllan N, Canfrán-Duque A, Aranda JF, Ramírez CM, Araldi E, 83. Zhang S, Hulver MW, McMillan RP, Cline MA, Gilbert ER. The pivotal et al. MicroRNA-148a regulates LDL receptor and ABCA1 expression to role of piruvate dehydrogenase kinases in metabolic flexibility. Nutr control circulating lipoprotein levels. Nat Med. 2015;21:1280–9. Metab. 2014;11:10. 60. Gastebois C, Chanon S, Rome S, Durand C, Pelascini E, Jalabert A, et al. 84. Piao L, Sidhu VK, Fang Y-H, Ryan JJ, Parikh KS, Hong Z, et al. FOXO1- Transition from physical activity to inactivity increases skeletal muscle miR- mediated upregulation of pyruvate dehydrogenase kinase-4 (PDK4) 148b content and triggers insulin resistance. Phys Rep. 2016;4:e12902. decreases glucose oxidation and impairs right ventricular function in 61. Rotllan N, Price N, Pati P, Goedeke L, Fernández-Hernando C. microRNAs in pulmonary hypertension: therapeutic benefits of dichloroacetate. J Mol lipoprotein metabolism and cardiometabolic disorders. Atherosclerosis. Med. 2013;91:333–46. 2016;246:352–60. 85. Gopal K, Saleme B, Al Batran R, Aburasayn H, Eshreif A, Ho KL, et al. FoxO1 62. Xu Q, Jiang Y, Yin Y, Li Q, He J, Jing Y, et al. A regulatory circuit of miR- regulates myocardial glucose oxidation rates via transcriptional control of 148a/152 and DNMT1 in modulating cell transformation and tumor pyruvate dehydrogenase kinase 4 expression. Am J Physiol Circ Physiol. angiogenesis through IGF-IR and IRS1. J Mol Cell Biol. 2013;5:3–13. 2017;313:H479–90. 63. Sarkar S, Dey BK, Dutta A. MiR-322/424 and −503 are induced during muscle 86. Ouyang W, Liao W, Luo CT, Yin N, Huse M, Kim MV, et al. Novel Foxo1- differentiation and promote cell cycle quiescence and differentiation by down- dependent transcriptional programs control Treg cell function. Nature. 2012; regulation of Cdc25A. Mol Biol Cell. 2010;21:2138–49. 491:554–9. 64. Jia H, Zhao Y, Li T, Zhang Y, Zhu D. miR-30e is negatively regulated by 87. Kim EH, Gasper DJ, Lee SH, Plisch EH, Svaren J, Suresh M. Bach2 regulates myostatin in skeletal muscle and is functionally related to fiber-type homeostasis of Foxp3+ regulatory T cells and protects against fatal lung composition. Acta Biochim Biophys Sin. 2017;49:392–9. disease in mice. J Immunol. 2014;192:985–95. 65. Zaragosi L-E, Wdziekonski B, Le Brigand K, Villageois P, Mari B, Waldmann R, 88. Itoh-Nakadai A, Matsumoto M, Kato H, Sasaki J, Uehara Y, Sato Y, et al. A et al. Small RNA sequencing reveals miR-642a-3p as a novel adipocyte- Bach2-Cebp gene regulatory network for the commitment of multipotent specific microRNA and miR-30 as a key regulator of human adipogenesis. hematopoietic progenitors. Cell Rep. 2017;18:2401–14. Genome Biol. 2011;12:R64. 89. Reznick J, Preston E, Wilks DL, Beale SM, Turner N, Cooney GJ. Altered 66. Shi C, Zhang M, Tong M, Yang L, Pang L, Chen L, et al. miR-148a is feeding differentially regulates circadian rhythms and energy associated with obesity and modulates adipocyte differentiation of metabolism in liver and muscle of rats. Biochim Biophys Acta Mol basis mesenchymal stem cells through Wnt signaling. Sci Rep. 2015;5:9930. Dis. 1832;2013:228–38. 67. Schweisgut J, Schutt C, Wüst S, Wietelmann A, Ghesquière B, Carmeliet P, 90. Dyar KA, Ciciliot S, Wright LE, Biensø RS, Tagliazucchi GM, Patel VR, et al. et al. Sex-specific, reciprocal regulation of ERα and miR-22 controls muscle Muscle insulin sensitivity and glucose metabolism are controlled by the lipid metabolism in male mice. EMBO J. 2017;36:1199–214. intrinsic muscle clock. Mol Metab. 2014;3:29–41. 68. Chen B, Tang H, Liu X, Liu P, Yang L, Xie X, et al. miR-22 as a prognostic 91. Yamaguchi S, Moseley AC, Almeda-Valdes P, Stromsdorfer KL, Franczyk MP, factor targets glucose transporter protein type 1 in breast cancer. Cancer Okunade AL, et al. Diurnal variation in PDK4 expression is associated with Lett. 2015;356:410–7. plasma free fatty acid availability in people. J Clin Endocrinol Metab. 2018; 69. Xiao D, Zhou T, Fu Y, Wang R, Zhang H, Li M, et al. MicroRNA-17 impairs 103:1068–76. glucose metabolism in insulin-resistant skeletal muscle via repressing 92. Everett LJ, Lazar MA. Nuclear receptor Rev-erbα: up, down, and all around. glucose transporter 4 expression. Eur J Pharmacol. 2018;838:170–6. Trends Endocrinol Metab. 2014;25:586–92. 70. Hudson NJ, Dalrymple BP, Reverter A. Beyond differential expression: the quest 93. Zhang X, Guo J, Wei X, Niu C, Jia M, Li Q, et al. Bach1: function, regulation, for causal mutations and effector molecules. BMC Genomics. 2012;13:356. and involvement in disease. Oxid Med Cell Longev. 2018;2018:1347969. 71. Óvilo C, Benítez R, Fernández A, Núñez Y, Ayuso M, Fernández AI, et al. 94. Li K, Qiu C, Sun P, Liu D, Wu T, Wang K, et al. Ets1-mediated acetylation of Longissimus dorsi transcriptome analysis of purebred and crossbred Iberian FoxO1 is critical for gluconeogenesis regulation during feed-fast cycles. Cell pigs differing in muscle characteristics. BMC Genomics. 2014;15:413. Rep. 2019;26:2998–3010.e5. 72. Li J, Yue Z, Xiong W, Sun P, You K, Wang J. TXNIP overexpression 95. Besse-Patin A, Jeromson S, Levesque-Damphousse P, Secco B, Laplante M, suppresses proliferation and induces apoptosis in SMMC7221 cells through Estall JL. PGC1A regulates the IRS1:IRS2 ratio during fasting to influence ROS generation and MAPK pathway activation. Oncol Rep. 2017;37:3369–76. hepatic metabolism downstream of insulin. Proc Natl Acad Sci U S A. 2019; 73. Ling L, Zhang S-H, Zhi L-D, Li H, Wen Q-K, Li G, et al. MicroRNA-30e 116:4285–90. promotes hepatocyte proliferation and inhibits apoptosis in cecal ligation 96. Laurie SM, Cain CC, Lienhard GE, Castle JD. The glucose transporter GluT4 and puncture-induced sepsis through the JAK/STAT signaling pathway by and secretory carrier membrane proteins (SCAMPs) colocalize in rat binding to FOSL2. Biomed Pharmacother. 2018;104:411–9. adipocytes and partially segregate during insulin stimulation. J Biol Chem. 74. Wrann CD, Eguchi J, Bozec A, Xu Z, Mikkelsen T, Gimble J, et al. FOSL2 1993;268:19110–7. promotes leptin gene expression in human and mouse adipocytes. J Clin 97. Yoshizumi S, Suzuki S, Hirai M, Hinokio Y, Yamada T, Yamada T, et al. Invest. 2012;122:1010–21. Increased hepatic expression of ganglioside-specific sialidase, NEU3, Mármol-Sánchez et al. Journal of Animal Science and Biotechnology (2020) 11:10 Page 21 of 21 improves insulin sensitivity and glucose tolerance in mice. Metabolism. 2007;56:420–9. 98. Jia Z, Moulson CL, Pei Z, Miner JH, Watkins PA. Fatty acid transport protein 4 is the principal very long chain fatty acyl-CoA synthetase in skin fibroblasts. J Biol Chem. 2007;282:20573–83. 99. Larkin JR, Zhang F, Godfrey L, Molostvov G, Zehnder D, Rabbani N, et al. Glucose-induced down regulation of thiamine transporters in the kidney proximal tubular epithelium produces thiamine insufficiency in diabetes. PLoS One. 2012;7:e53175. 100. Rey G, Valekunja UK, Feeney KA, Wulund L, Milev NB, Stangherlin A, et al. The pentose phosphate pathway regulates the circadian clock. Cell Metab. 2016;24:462–73. 101. Baenziger JE, Hénault CM, Therien JPD, Sun J. Nicotinic acetylcholine receptor–lipid interactions: mechanistic insight and biological function. Biochim Biophys Acta Biomembr. 1848;2015:1806–17. 102. Liu LB, Chen XD, Zhou XY, Zhu Q. The Wnt antagonist and secreted frizzled- related protein 5: Implications on lipid metabolism, inflammation, and type 2 diabetes mellitus. Biosci Rep. 2018;38:BSR20180011. 103. Casas E, Shackelford SD, Keele JW, Stone RT, Kappes SM, Koohmaraie M. Quantitative trait loci affecting growth and carcass composition of cattle segregating alternate forms of myostatin. J Anim Sci. 2000;78:560. 104. Nalaila SM, Stothard P, Moore SS, Li C, Wang Z. Whole-genome QTL scan for ultrasound and carcass merit traits in beef cattle using Bayesian shrinkage method. J Anim Breed Genet. 2012;129:107–19. 105. Schering L, Albrecht E, Komolka K, Kühn C, Maak S. Increased expression of thyroid hormone responsive protein (THRSP) is the result but not the cause of higher intramuscular fat content in cattle. Int J Biol Sci. 2017;13:532–44. 106. Ortega-Prieto P, Postic C. Carbohydrate sensing through the transcription factor ChREBP. Front Genet. 2019;10:472.