Correlated mRNAs and miRNAs from co-expression and regulatory networks affect porcine muscle and finally meat properties

Meat quality and carcass traits phenotypes

Gene co-expression network construction for mRNA

Module-trait associations of mRNA

Sets of genes (modules) with common expression patterns that were associated with particular traits were identified based on the correlation between ME and organismal phenotype. We identified five modules that significantly associated with meat quality. Modules dark-turquoise and orange were correlated positively to pH traits and negatively to drip loss (ME_{[dark-turquoise]}: pH24MLD r = 0.34, p = 5.3 × 10^-7, DL r = -0.19, p = 5.6 × 10^-3; ME_[orange]: pH24MLD r = 0.32, p = 3.7 × 10^-6, DL r = -0.31, p = 5.8 ×10^-6) (Figure 1). Module dark-turquoise (31 annotated genes) was highly enriched for genes belonging to the cluster “glucose metabolic process” (GO: 0006006) and the KEGG-pathway “insulin signaling” with an enrichment score (ES) of 2.65. Module orange (26 annotated genes) was enriched for transcripts of the functional annotation clusters “response to wounding”, “defense response” and “inflammatory response” (ES = 2.42). Modules red, black, and tan were correlated negatively to pH traits and positively to drip loss (ME_[red]: pH45MLD r = -0.22, p = 1.8 × 10^-3, DL r = 0.20, p = 3.9 × 10^-3; ME_[black]: pH45MLD r = -.23, p = 8.8 × 10^-4, DL r = 0.16, p = 1.8 × 10^-2; ME_[tan]: pH45MLD r = -0.22, p = 1.4 × 10^-3, DL r = 0.19, p = 6.9 × 10^-3) (Figure 1). Further, modules red (315 annotated genes), black (436 annotated genes), and tan (154 annotated genes) were enriched for genes of the top functional annotation clusters of “mitochondrial ribosome”, “mitochondrion”, and “extracellular matrix” with ES of 10.23, 15.15, and 27.05, respectively. Only one module (ME_{[dark-orange]}) showed association with traits related to fatness (Figure 1).

Co-expression networks and module-trait associations for miRNA

Transcriptional networks of muscle microRNAs were studied with a WGCNA using miRNA expression data on M. longissimus dorsi from the same animals as above. Residuals, derived from the analysis of expression levels after correction for systematic effects according to the mixed model, were used for constructing the muscle miRNA transcriptional networks, i.e., modules. Nine modules were identified (Figure 2). Only 2 modules were associated with meat quality at a significance level of p < 0.05. Module purple was correlated positively to LF24MLD at p = 0.03 and negatively to pH24MLD and IMP24MLD at p = 0.04 and p = 0.03, respectively. Module purple consisted of 8 miRNA families (miR-17, miR-30. miR221, miR-185, miR-324, miR362, miR-500, and miR-542). Module blue was positively correlated to pH45MLD and IMP24MLD at p = 0.02 and p = 0.01, respectively. Module blue comprised 29 miRNA families (let-7, miR-15, miR-17, miR-31, miR-95, miR-103, miR-105, miR-122, miR-124, miR-130, miR-138, miR-154, miR-184, miR-185, miR-197, miR-202, miR-204, miR-212, miR-214, miR-320, miR-326, miR-335, miR-346, miR-383, miR-467, miR-491, miR-744, miR-1224, and miR-1296).

Individual miRNA expression profiles and correlated traits

Endogenous correlation of expression profiles between miRNA and mRNA

Integration of negative correlation of miRNA and mRNA with module-trait association and target prediction

A total of 1,928 pairs of miRNAs and mRNAs that showed negative correlations at p-values ≤ 8.97 × 10^-5 (FDR = 0.1) belonged to 1,073 mRNA probe sets (929 gene) and 264 miRNA probe sets (96 miRNA families). Of these, 286 pairs were assigned to modules dark-turquoise, red, black, and tan, which showed correlation with traits related to meat quality (Figure 3). However, no genes in module orange were negatively correlated with miRNA at FDR < 0.1. Only one gene (CREM) in the module dark-turquoise was negatively correlated with miR-153 (r = -0.34 p = 2.11x10^-06 FDR = 0.02). In module black, 69 genes were negatively correlated with 21 miRNA families, totaling 140 miRNA-mRNA pairs. In total, 101 and 43 pairs of miRNAs and mRNAs were identified in modules red (22 miRNAs and 52 genes) and tan (17 miRNAs and 26 genes). Known genes belonging to modules dark-turquoise, red, black, and tan and their negatively-correlated miRNAs are shown in Table 5. Out of 1,928 pairs of miRNAs and mRNAs that showed negative correlations at p-values ≤ 8.97 × 10^-5 (FDR = 0.1), 62 pairs were assigned to modules blue and purple, which correlated with meat quality (Figure 4). These 62 pairs of miRNAs and mRNAs belonged to 14 miRNA families and 40 genes (Table 6).

Further, TargetScan and RNAhybrid were used to scan miRNA and mRNA sequences (porcine RefSeq), to obtain additional evidence for their functional link; these sequences corresponded to 1,928 pairs of negatively-correlated miRNAs and mRNAs. In total, 474 pairs of miRNA and mRNA were confirmed by either of the two in silico prediction methods: 331, 195, and 32 miRNA-mRNA pairs were predicted by RNAhybrid, TargetScan, or both methods, respectively. The 474 miRNA-mRNA pairs covered 121 probe sets of miRNAs (65 miRNA families) and 331 targets probe sets (297 genes). When focusing on genes of the trait-correlated modules dark-turquoise, orange, red, black, and tan, 73 out 286 pairs of miRNAs and mRNAs were confirmed with at least one in silico method. These 73 pairs comprised 26 miRNA families and 51 genes (Tables 5 and 6 in bold).

mRNA abundance and co-expression networks linked to muscle and meat properties

We used transcriptional network analysis to identify co-expression modules (dark-turquoise, orange, red, black, and tan) that correlated to meat quality phenotype. These modules were highly enriched for genes involved in “extracellular matrix”, “glucose metabolic process”, and “mitochondrion” (“oxidative phosphorylation” KEGG_PATHWAY); i.e. processes affecting structural and metabolic properties.

A dominant role of mitochondria is the production of ATP by oxidative phosphorylation that depends on oxygen supply. When oxygen is limited (post mortem or during prolonged vigorous exercise) the glucose metabolism occurs by anaerobic respiration, a process that is independent of the mitochondria. A shift from aerobic to anaerobic metabolism - favouring the production of lactic acid - results in a pH decline post mortem and thereby influence the meat quality [29]. So the biological process of mitochondria as well as the way of glucose metabolism play a significant role in the muscle cell and finally impact on meat quality. Indeed, mitochondrial dysfunction resulting in decreased cellular energy production is also responsible for a variety of human myopathies and cardiomyopathies [30–33].

Genes assigned to the GO category of “extracellular matrix” encode proteins belonging to the myofibrillar scaffold. The characteristics of the myofibrillar scaffold and the kinetics of their fragmentation were associated with tenderness and water-holding capacity of meat. In particular, the proteolysis of muscle proteins affects the shrinkage of myofibrils, the development of pores in the cell membranes, so called drip channels, and the non-covalent binding of water molecules [34, 35]. Collagens are major constituents of the extracellular matrix (ECM). In our study there were many collagen genes that are reported to be correlated with various muscle disorders [36–39]. For example, collagen type VI (COLVI), an important component of skeletal muscle ECM, is involved in maintaining tissue integrity [40]. Col6a1^-/- mice show a complete absence of collagen VI chains and display a myopathic phenotype, abnormal mitochondria, and increased apoptosis of muscle fibres [41, 42].

miRNA abundance and co-expression networks linked to muscle and meat properties

Our finding of the relevance of mitochondrial metabolic pathways, including oxidative phosphorylation, and muscle structural protein composition to post mortem processes affecting meat quality is in line with our previous results obtained in other populations [1, 6, 7]. Additionally, while the previous studies focused on trait-associated mRNA expression, here another hierarchical level in the regulatory network relevant to processes occurring under conditions of insufficient oxygen, energy, and nutrient supplies is provided. In fact, miRNA was integrated into this study as a regulator molecule of muscle transcripts. miRNAs with identical seed sequences (the same family) [43] or that are closely located on the same chromosome (the same cluster) [44, 45] have similar expression trends. This was confirmed in our study, where most of the modules consisted of the same families of miRNAs or miRNAs located on the same chromosome. In this study, marginal association of miRNA co-expression modules to organismal traits was found compared to mRNA co-expression modules. On the one hand, co-expressed miRNAs of the same family or cluster might not regulate the same trait. On the other hand, this may be caused by indirect regulation of organismal traits by miRNAs via their effect on mRNA transcripts. Accordingly, individual miRNA correlation to phenotypes was also considered. Recent studies have revealed key roles for miRNAs in the regulation of skeletal muscle differentiation, and changes in miRNA expression are associated with various skeletal muscle disorders [46–48]. In this study, several miRNAs were correlated with carcass and meat quality traits. This includes miR-221, previously identified in studies of myotube maturation and in the maintenance of the myofibrillar organization [49] and found to contribute to muscle pathogenetic mechanisms [50]. Interestingly, miR-133, which showed highest correlation with pH24MLD and pH24MSM, has been widely studied for roles in the regulation of skeletal muscle development, including in proliferation and myogenesis [10, 51] as well as muscle disorders [47]. Recently, a study reported that mice with genetic deletions of miR-133a-1 and miR-133a-2 developed adult-onset centronuclear myopathy in type II (fast-twitch) myofibres, which was accompanied by impaired mitochondrial function, fast-to-slow myofibre conversion, and disarray of muscle triads [52]. These are changes of muscle structure and metabolism that also impact meat quality. In addition to its well established role in translation (Table 4), miR-23a was also recently identified as a key regulator of skeletal muscle differentiation and is predicted to target multiple adult fast myosin heavy chain (Myh) genes, including Myh1, 2, and 4[53]. For fat traits, miR-103 and miR-107 were highly correlated (Table 3). This is consistent with previous reports of miR-103 being involved in adipogenesis, lipid metabolism, and adipocyte differentiation [54, 55] and of miR-103/107 being involved in glucose homeostasis and insulin sensitivity [56].

In this study we showed for the first time more complex correlations among miRNAs and between miRNAs and post mortem organismal phenotypes in swine, while also confirming previous studies in human and mouse muscle as well as C2C12 myoblasts.

Links between miRNA and mRNA that relate to muscle and meat properties

This study also sought to evaluate to what extent the co-expression modules of trait-associated mRNAs are themselves regulated by miRNAs. A regulatory link between miRNA and mRNA and a functional link to the organismal phenotype was suspected if (1) the mRNA belonged to either one of the co-expression modules associated with the traits (i.e., ME_{dark-turquoise}, ME_orange, ME_red, and ME_tan), (2) mRNA abundance was significantly negatively correlated with its miRNA regulator, and (3) the mRNA was predicted to be a target gene of the respective miRNA. Therefore, RefSeqs of the genes with present calls from the 3′-IVT-Affymetrix arrays were explored to predict the targets of miRNAs by either seed sequence complementarity [57] or by thermodynamics-based modeling of RNA:RNA duplex interactions [58]. Currently, no publicly-accessible database covers porcine miRNAs and their predicted target genes. Moreover, annotation of porcine genes is not yet finalized. Accordingly, the target predictions should be interpreted with caution.

Interestingly, no miRNA regulator was identified by negative correlation or target prediction in the module orange, which was enriched for genes related to “response to wounding”, “inflammatory response”, and “defense response”. Genes assigned to biofunctions related to response to exogenes stimuli, change their transcription rate immediately due to many factors. As suggested by finding no correlated miRNAs – their regulation of transcription may occur without major involvement of miRNA. However, many genes in the module orange were previously confirmed as transcriptional regulators in myogenesis or were located in QTL regions for muscle fiber traits like BTG2, EGR1, ANKRDS1 and FOS [59, 60]. Interestingly, the genes in module orange like Egr1, FOS and JUN that are associated with oxidative stress response were found upregulated in muscle in response to mechanical ventilation and immobilization in a porcine model for critical illness myopathy (CIM) [61, 62].

Co-expressed genes in module dark-turquoise were significantly associated with meat quality and based on the current knowledge of gene functions some links among them are suggested to be relevant For example, one member of this module, CREM, is a transcription factor binding to cAMP-responsive elements (CREs) in the promoters of various genes. This transcription factor plays important roles in various organismal functions [63–67]. Crem inactivation or knockout has been shown to increase the rate of apoptosis in testis tissue [68, 69]. The main cellular change associated with apoptosis processes also occur during post-mortem [70]. Post mortem biochemical processes in muscle lead to pH decline. A high expression of CREM being positively correlated with pH at 24 hours, may indicate a slowdown of apoptosis related post mortem processes paralleling anaerobic metabolic processes that led to a decrease of pH. This indicates that the abundance of CREM transcripts in muscle plays a significant role in meat quality. Further, CREM and miR-153 were highly negatively correlated (FDR < 0.01) which is known to induce apoptosis in a glioblastoma cell line DBTRG-05MG [71]. Thus a consistent link of effects on apoptosis and mRNA-miRNA interaction is obvious. miR-153 also inhibits the protein kinase B (PKB/Akt) pathway by reducing the protein level of insulin receptor substrate-2 (IRS2) [72]. As recently shown, miR-135a targets IRS2 levels by binding to its 3′UTR and this interaction regulates skeletal muscle insulin signaling [73]. Insulin signaling plays a pivotal role in the regulation of glucose uptake by skeletal muscle [74]. The glucose uptake in skeletal muscle has large effects on meat characteristics [75, 76]. In our study, Irs-2 also belonged to module dark-turquoise and was negatively correlated with waterholding capacity related traits like DL (p = 0.004) and positively with pH (pH24MLD, p = 7.05E-07; pH24MSM, p = 1.30E-03; pH45MLD, p = 3.39E-03). Thus another plausible functional link of members of the module dark-turquoise and miRNAs can be shown.

Mitochondria supply energy for physiological functions and play a significant role in the regulation of other cellular events including apoptosis, calcium homeostasis, and production of reactive oxygen species. Mitochondrial metabolism is affected by miRNA regulation [77]. Here, we found many miRNAs being negatively correlated to target genes of modules red, black, which were enriched for genes related to mitochondrial pathways. Indeed, miR-338 modulates energy metabolism, oxidative phosphorylation, and mitochondrial functions [78, 79], and miR-15b, -16, -195, and -424 decrease cellular ATP levels in cardiomyocytes [80]. Additionally, miR-181c can enter and target the mitochondrial genome, ultimately causing electron transport chain complex IV remodeling and mitochondrial dysfunction [81]. Here we found miR-181 was correlated with DL and pH24MLD as well as S100 calcium binding protein A6 (S100A6). Further, mitochondrial genes like Kif1b, Atp6v1g1, Atp5c1, Park7, Chchd2, Ruvbl1, Mrps7, and Mrpl48 were highly negatively correlated with, and some of them were predicted as targets of, miR-15, -154, -17, -503, -214, -330, -342, and -375.

Many genes of modul tan were assigned in the GO category “extracellular matrix” including Col1a1, Col1a2, Col3a1, Col5a1, Col5a2, Col6a1, Col6a3, Col12a1, Col14a1, Crispld2, Ddn, Fbn1, Nid1, Nid2, Ogn, and Postn. These genes were negatively correlated with, and some of them were predicted as targets of, miRNA. Interestingly, Col6a1 was found as a target for many miRNAs including miR-132, miR-205, miR-320, miR-330, miR-375, miR-423, and miR-425. Moreover, miR-205 was identified as a regulator of Col1a1, Maged1, and Dpysl3. Gandellini et al. (2012) [82] reported that miR-205 controls the deposition of laminin-332 and its receptor integrin-β4 as well as participates in a network involving ΔNp63α, which is essential for maintenance of the basement membrane in prostate epithelium. Similarly, other miRNAs, including miR-29, miR-133, and miR-30, are involved in the regulation of development and maintenance of extracellular matrix of bone and muscle [83, 84].

Much evidence suggests that a group of miRNAs (cluster and family) may contribute to the regulation of a set of common targets [46, 85, 86], and are, therefore, associated with phenotypes. WGCNA was used here to group miRNA products and revealed 2 modules associated with meat quality (purple and blue). The miRNA from these modules were negatively correlated to mRNAs, and some of these were predicted as targets. Most miRNAs in module purple were related to genes in the categories “ubiquitin” or “protein catabolic process” (Ube2e1, Ube2b1, Ubap1, Igfbp5, Herc1, Psmb3, Flj25076). Differential expression of genes of the ubiquitin system depending on muscle and meat quality was previously shown; only recently also association of genes of the ubiquitin system with meat quality was reported [1, 87]. In particular, miR-324 was highly negatively correlated and predicted to target Ube2e1, Ube2v1, Bap1, Lpf and Cic.

Most previous muscle mRNA and miRNA expression studies focused on cardiac muscle or skeletal muscle injury [50, 88]. In these injuries, the degree of damage results from an imbalance of energy, nutrients, and oxygen supply-and-demand in muscle cells. Similarly, nutrient, energy, and oxygen depletion occurs post mortem. Many changes in expression associated with muscle injury would therefore overlap with post mortem processes in conversion of muscle to meat, and vice versa. In this regard, functional annotation of mRNA co-expression and trait-correlated expression identified key post mortem pathways and functions, including glucose metabolic process, mitochondrial metabolic pathways, and muscle structural components, involved in muscle-to-meat conversion that will be relevant to muscle injury as well.

Animals, tissue collection, and phenotyping

Animal care and tissue collection procedures followed the guidelines of the German Law of Animal Protection, and the experimental protocol was approved by the Animal Care Committee of the FBN. This study was based on trait measurement and expression profile association analyses done with 207 performance-tested pigs from commercial herds of the crossbreed Pietrain × (German Large White × German Landrace). Animals were raised and slaughtered under standardized conditions in the experimental facilities of the Leibniz Institute for Farm Animal Biology (FBN). Sample collection was performed thoroughly after exsanguination, tissue samples were rapidly dissected, snap-frozen in liquid nitrogen and stored at -80°C. The average age of the pigs at sampling was ~180 days. Technological parameters of meat quality, i.e., pH-value, conductivity, and colour, were measured by using Star-series equipment (Rudolf Matthaeus Company, Germany). Measures of pH and conductivity were at 45 min post mortem (pH45) and 24 h post mortem (pH24), in both M. longissimus dorsi between 13th/14th rib (pH45MLD, pH24MLD, LF45MLD, LF24MLD) and the ham (M. semimembranosus) (symbol: pH24MSM, LF24MSM). Muscle colour was measured at 24 h post mortem by Opto-Star (Matthaeus, Klausa, Germany). Drip loss was scored based on a bag method with a size-standardized sample from the M. longissimus dorsi collected at 24 h post mortem and weighed, suspended in a plastic bag, held at 4°C for 48 h, and re-weighed [91]. To determine cooking loss, a loin cube was taken from the M. longissimus dorsi, weighed, placed in a polyethylene bag, and incubated in water at 75°C for 50 minutes. The bag was then immersed in flowing water at room temperature for 30 minutes, and the solid portion was re-weighed. Thawing loss was determined similarly after at least 24 h freezing at -20°C. Drip loss, cooking loss, and thawing loss were calculated as a percentage of weight loss based on the start weight of a sample. Shear force was measured using Instron-4310 equipment, and average values of four replicates were used for analyses.

Customized miRNA microarrays design

Our custom porcine miRNA array was designed from 284 pig miRNAs obtained from the miRBASE (miRBase 14.0). Because miRNAs are highly conserved between closely-related species [10], we could predict novel porcine miRNA candidates by inter-species alignments requiring 100% mature miRNA similarity. Accordingly, we used previously known miRNA sequences from humans and mice, as well as other species, to perform BLAST searches against the porcine genome database porcine; 391 miRNA candidates were identified. In total, 675 miRNAs probe sets, corresponding to 513 unique sequences belonging to 159 miRNA families, were used for hybridisation with the target samples described above. Microarray data related to all samples were deposited in the Gene Expression Omnibus public repository (GEO accession number: GSE41294: GSM1013731-GSM1013920).

Customized microarrays, pre-processing, and normalization of miRNA

Total miRNA was isolated with Qiagen miReasy Mini kit and RNeasy MinElute Clean up kit (Qiagen, Hilden, Germany) according to manufacturer’s protocol for small RNA. Quality and quantity of isolated total RNA and miRNA were determined using an Agilent 2100 Bioanalyzer for RNA (Agilent Technologies, Santa Clara, CA). Affymetrix customized microarrays from our porcine miRNA candidate dataset were used. Targets for hybridisation were prepared from miRNA with the FlashTag™ Biotin RNA Labeling Kit for Affymetrix GeneChip miRNA arrays (Genisphere, Hatfield, PA, USA) according to manufacturer’s recommendations. Briefly, 250 ng of miRNA of each individual were poly(A)-tailed using ATP–poly-A-Polymerase, then FlashTag Biotin end-labelled. After hybridisation of biotin-labelled complementary RNA, chips were washed and processed to detect biotin-containing transcripts by Streptavidin-PE (Phycoerythrin) conjugate, then were scanned on GeneChip scanner 3000 7G (Affymetrix, Santa Clara, US). Data were extracted from the images, and spots were quantified and processed by quality filtering. Hybridisation quality was assessed in all samples by using JMP Genomics 5 utilising Robust Multi-array Average (RMA) background correction and log2 transformations. To acquire the expression value, data were normalized between chips using the quantile normalization method.

Whole-genome expression profiling (mRNA)

Gene expression profiling of M. longissimus dorsi samples of pigs was conducted with the same animals (207) as for miRNA. In brief, total RNA was isolated using TRI Reagent (Sigma, Taufkirchen, Germany) and used for target preparation for microarray hybridisation. According to Affymetrix protocols, 500 ng of total RNA were reverse-transcribed into cDNA, transcribed into cRNA, and labelled using Affymetrix One cycle synthesis and labelling kit (Affymetrix, UK) to prepare antisense biotinylated RNA targets. Quality of hybridisation was assessed in all samples following manufacturer’s recommendations. Data were analysed with the Affymetrix GCOS 1.1.1 software, using global scaling to a target signal of 500. Data were processed with MAS5.0 to generate cell intensity files (present or absent). Quantitative expression levels of present transcripts were estimated using the PLIER algorithm (Probe Logarithmic Intensity Error) for normalization that was implemented in Expression Console (Affymetrix). Based on BLAST comparison of the Affymetrix porcine target sequences with the porcine genome sequence (Ensembl_Sscrofa_10), 20,689 of the 24,123 probe sets on the Affymetrix Porcine GeneChip were localized and annotated [92, 93]. Microarray data related to all samples were deposited in the Gene Expression Omnibus public repository [GEO accession number: GSE32112: GSM796045-GSM796251].

Pre-processing of phenotype and expression data

Phenotypes and expression levels were adjusted for systematic effects by analysis of variance performed with the procedure “Mixed” of the SAS software package (SAS version 9.1 SAS Institute, Cary, NC) before analysing their correlation and by using co-expression network. Sex and RYR genotype was used as a fixed effect, “sire” and “slaughter day” as random effects, and “carcass weight” as a covariate. Subsequently, the residuals of log2-transformed expression intensities (miRNA and mRNA) and muscle phenotype were used for further analysis.

Pearson correlation of miRNA expression level and gene expression level was calculated using 190 individuals; correction for multiple testing was done by controlling the FDR level (q-value according to Storey and Tibshirani, 2003, [94]) at 10%.

Weighted gene Co-expression network analysis (WGCNA)

A weighted gene co-expression network was constructed for 207 muscle biopsies using the blockwise Modules function from the WGCNA package in R [23]. Residuals of gene expression, after correcting the effect, were used for WGCNA. The blockwise Modules function allows the entire dataset of 11,191 probe sets by mRNA and 675 miRNA to be utilised in the construction of the weighted gene co-expression network.

Extremely outlying individuals were removed from the following analysis based on hierarchically clustered using the average linkage function, and common Euclidean distance. Pearson correlation matrix of all gene-gene comparisons were calculated across all microarrays. Adjacency matrix was then calculated using the correlation matrix of the expression sets. Finally, the Topological Overlap Matrix (TOM) was converted from the adjacency matrix and used to derive a TOM-based distance matrix for the hierarchical clustering of expressions. In the next step, modules of expression profiles, (i.e. sets of genes with high topological overlap) were formed based on hierarchical clustering, with empirically specified minimal module size (30 for gene expression, 10 for miRNA).

where x _i and x _j are the gene expression profiles of the i-th and j-th gene and α _ij is the adjacency. TOM based distance matrix is a robust and powerful measurement in building co-expression network. Selection of appropriate value for the power β were derived according to the pickSoftThresholding function of the WGCNA package [23, 95]. Accordingly, by manually inspecting the fit of the scale free topology model with the candidate β-values for each set of expression profile, minimal β-values giving a coefficient of determination R² higher than 90% were adopted.

Modules were further merged based on the dissimilarity between their “eigengenes”, which were defined as the first principal component of each module. Genes that were not assigned to another module were assigned to the grey module. A threshold of 0.2 for the dissimilarity as recommended by the WGCNA author was used. Module–trait associations were estimated using the correlation between the module eigengene and the phenotype, which allows easy identification of expression set (module) highly correlated to the phenotype. For each expression profile, Gene Significance (GS) was calculated as the absolute value of the correlation between expression profile and each trait; module membership (MM) was defined as the correlation of expression profile and each module eigengene, enabling further identification of key players in the regulation network.

Gene ontology and pathway enrichment analysis

We performed a gene ontology (GO) enrichment analysis for network modules using the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/[96, 97]. Each analysis was performed using the functional annotation clustering option. Functional annotation clustering combines single categories with a significant overlap in gene content and assigns an enrichment score (ES, defined as the -log10 of the geometric mean of unadjusted p-value s for each single term in the cluster) to each cluster, making interpretation of the results more straightforward. To assess the significance of functional clusters, we created 22 sets of 11,191 probe sets corresponding to 8,036 genes (size of the average module identified in this study).

Predicting porcine targets of miRNAs by RNA hybrid and TargetScan

We used two methods to predict the targets of porcine miRNA. First, we predicted targets using the computational software RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid), which detects the most energetically favourable hybridisation sites of a small RNA within a large RNA [98]. Here, we tested the miRNA probe sets with the following parameters: number of hits per target = 1, energy cutoff = -25 kcal/mol, and maximal internal or bulge loop size per side = 4. Most targets found were located on the 3′-UTR of genes. Second, TargetScan (http://www.targetscan.org) was used to detect target gene candidates based on seed complementarity on UTR database 6.0 and our porcine RefSeq transcript with our miRNA seed sequence [99]. TargetScan was applied considering both conserved and non-conserved targets. The porcine RefSeq transcripts, which derived from 11,191 probe sets that showed consistent expression in porcine muscle, were used as input targets for RNAhybrid and TargetScan.