BACKGROUND
micrornas  are endogenous small non-coding rnas ~22 nt in length that play critical roles in diverse biological processes via epigenetic regulation of gene expression.
precursor mirnas  are initially generated in nucleus and processed into an approximately 70 nt long stem-loop structure.
it is then exported to cytoplasm and processed by dicer to generate mirna/mirna duplexes.
one strand of which is incorporated with agonaute to form rna-induced silencing complex  that targets mrnas via base-pair complementary, typically to their 3′ untranslated regions  or cds and downregulates gene expression by either degradation of mrna or repression of translation, while the other strand is usually discarded  <cit> .
skeletal muscle is highly metabolically active and valuable for meat-producing animals.
slow-twitch-oxidative , fast-twitch-oxidative , and fast-twitch-glycolytic  fiber were the three major muscle fiber types in pigs.
muscle fibers have strong association with muscle metabolic activities and meat quality such as tenderness, juiciness and color.
muscle containing a high proportion of oxidative fibers is often associated with higher fat content, oxidative enzyme activities and mitochondrial density ; a high ratio of ftg fibers is associated with high glycolytic enzyme activities.
previous research has identified several mirnas are associated with meat quality such as mir- <dig>  mir- <dig> and mir- <dig> etc.
in porcine skeletal muscle  <cit> .
the polymorphisms in the porcine mir- <dig> and mir- <dig> are proposed as a genetic factor affecting muscle fibers and meat quality traits  <cit> .
since the critical roles of mirnas such as myogenesis, adipogenesis and muscle development have been discovered in pig skeletal muscle , the understanding of the mirna regulation in metabolic properties of skeletal muscle fibers could be key to improvement of meat quality  <cit> .
mir- <dig> and mir- <dig> could regulate the gene expression of oxidative phosphorylation  machinery including complex iv subunits cox <dig>  coxiv and atp synthase subunits atp5g <dig>  <cit> .
mir-15a and mir-15b modulate the cellular atp levels  <cit> .
mir- <dig> regulates the fatty acid oxidation and mitochondrial biogenesis through targeting peroxisome proliferator-activated receptor-gamma coactivator-1alpha   <cit> .
with those mirnas identified, the lack of a comprehensive and systematic mirna profiling associated with energy metabolism of skeletal muscle remains unraveled.
our previous research on muscle transcriptional profile has revealed numerous biological pathways significantly associated with muscle fiber type, mitochondrial respiratory activity and metabolic enzymes  <cit> .
it is of interesting to further investigate how mirnas are involved in energy metabolism by fine-turning gene expression.
in the present study, the mirna transcriptome profiling of longissimus muscle  samples obtained 24 h before slaughter of two pig breeds duroc and pietrain exhibiting divergent meat quality and muscle phenotypes provided a comprehensive insight into the discovery of mirnas associated with muscle fiber, mitochondrial respiratory activity and metabolic enzyme activity.
muscle of duroc pigs contains higher percentage of sto fibers, mitochondrial respiratory activity and higher fat content to improve the tenderness and juiciness of the meat.
in comparison, pinn pigs are more muscular and favorable for meat industry.
their skeletal muscles are leaner and contain more ftg fibers .
hence duroc and pietrain pigs are great models to study energy metabolism of skeletal muscle.
the mirna and mrna expression profile was then integrated based on their pairwise correlations and computational target prediction to construct the regulatory mirna-mrna networks which could potentially affect metabolic properties of skeletal muscle and hence meat quality.
the illumination of mirna-based regulatory metabolism could enrich our knowledge of the roles of mirnas in achieving phenotypic diversity of skeletal muscle in different breeds.
methods
animals and sample collection
the experiment and muscle biopsy collection were approved and authorized by the german and european animal welfare regulations for animal husbandry, transport, and slaughter .
all experimental procedures, including animal care and tissue sample collection, followed guidelines for safeguarding and good scientific practice in accordance with the german law of animal protection, officially authorized by the animal care committee and authorities .
as previously described , duroc and pietrain  pigs were raised until 180 days of age.
to rule out the effects of the maglinant hyperthermia syndrome  locus, only muscle samples from mhs-negative genotype pigs were investigated.
muscle biopsies were collected from five female and male pigs of each breed  for phenotypic measurements  .
biopsies were collected 24 h before slaughter from the longissimus muscle between 13th and 14th thoracic vertebrae.
muscle samples were frozen in liquid nitrogen and stored at −80 °c until analysis.
rna isolation
small rnas were isolated and enriched from longissimus muscle biopsies using a mireasy mini kit and an rneasy minelute cleanup kit  according to manufacturer’s protocols.
the quality and quantity of small rna were assessed with an agilent  <dig> bioanalyzer  using an agilent small rna kit.
microrna microarray analysis
the affymetrix genechip mirna  <dig>  array  was used to determine the mirna expression profile of the lm 24 h ante mortem of duroc and pinn pigs.
it is comprised of  <dig>  entries representing hairpin precursor, total probe set  <dig>  for detection of most mirna from  <dig> species .
200 ng of small rna were used in sample preparation using a flashtag biotin hsr rna labeling kit  for affymetrix genechip mirna arrays.
each of the labeled rna  was then hybridized for 16 h to an affymetrix genechip mirna array according to manufacturer’s recommendation , washed and stained in the affymetrix fluidics station  <dig>  and scanned on the affymetrix g <dig> gene array scanner.
expression console software was used for robust multichip average  normalization and the detection of present mirnas by applying the dabg  algorithm.
further filtering was done by excluding probes that were present in less than 70 % of the samples within each breed and annotated mirnas that had sequence greater than or equal to 30 nt in length.
three thousand five hundred eighty seven probes passed the quality filtering and were used for further analysis.
the availability of expression data are in the gene expression omnibus public repository with the geo accession number gse80198: gsm2120718-gsm <dig>
statistics and bioinformatics analyses
differential expression analysis for mirna was performed using the anova procedure in jmp genomics  <dig>  breed was treated as a fixed effect.
false discovery rate  was used to control an error rate of a multiple-hypothesis testing according to benjamin & hochberg  <cit> .
we used our previous microarray-based mrna expression data to integrate with the differentially expressed mirnas and scan for potential target genes.
pearson correlation of mirna and mrna expression levels was calculated.
both rnahybrid  <dig> .
<dig> and targetscan  <dig>  were used to predict targets of mirnas.
rnahybrid  is computational software that detects the most energetically favorable hybridization sites of a small rna within a large rna  <cit> .
the mirna probe sets were tested with the following parameters: number of hits per target =  <dig> and energy cutoff = −25 kcal/mol and maximal internal or bulge loop size per side =  <dig>  targetscan  was used to predict the target gene candidates based on complementarity of the mirna seed sequence  and target binding site on the 5′ utr, 3′ utr and protein coding region of the porcine mrna sequences   <cit> .
hafner et al.
and chi et al.
have found that argonaute-bound target sites in coding sequences  are as numerous as those located in 3′ utr in both hek <dig> cells   <cit>  and mouse brain   <cit> .
other research suggests that mirna target sites in 3′utr are more efficient at triggering mrnas degradation while cds and 5′utr located sites can effectively repress translation  <cit> .
xu et al.
develop novel computational approach for target prediction with sites located along the entire gene sequences to increase the percentage of true positive targets  <cit> .
hence 5′ utr, cds and 3′utr were included in this study to improve the sensitivity of mirna target identification and avoid a substantial number of missing targets.
transcripts that negatively correlated with mirna and predicted as potential targets were further passed to functional analysis.
ipa software  was used for functional analysis of potential mirna target genes.
it categorizes genes based on annotated gene functions and statistically tests for over-representation of functional terms within the gene list using fisher’s exact test.
the mirna-mrna regulatory networks were visualized using cytoscape v <dig> .
<dig>   <cit> .
quantitative real time pcr  for microrna microarray  validation
four mirnas  related to energy metabolism were validated by qpcr of each individual sample .
the same samples were used for both qpcr validation and mirna-chips.
the cdna was synthesized from 250 ng isolated mirnas using an ncodetm vilotm mirna cdna synthesis kit  according to manufacturer’s protocols.
real-time pcr was performed using the express sybr greenertm mirna qrt-pcr kit with premixed rox  according to manufacturer’s protocols.
all measurements were performed in duplicates.
the thermal parameters were 50 °c for 2 min, 95 °c for 2 min, followed by 45 cycles of 95 °c for 15 s and 60 °c for 1 min.
the universal qpcr primer was provided in the kit and the mirna-specific forward primers were designed for the mirnas of interest.
the designed primer sequence information is accessible in additional file  <dig>  geometric mean of the 5s and u <dig> transcription levels was used as an internal standard to normalize the mirna expression value.
correlation coefficient analysis between the michip and qpcr data was calculated using sas  <dig>  .
RESULTS
differentially expressed mirnas between duroc and pietrain
out of  <dig> probes quality-filtered probes,  <dig> probes belonged to  <dig> unique mature mirna sequences were differentially expressed  between duroc and pinn using anova on jmp genomics  <dig>  among these,  <dig> probes belonged to  <dig> mature mirna sequences were up-regulated in duroc pigs, while  <dig> probes belonged to  <dig> different mature mirna sequences were up-regulated in pinn pigs .
among these, mir- <dig>  mir- <dig> and mir- <dig> were the top three upregulated mirnas in duroc pigs with fold change ranging from  <dig>  to  <dig> .
whereas the top three upregulated mirnas in pinn were mir- <dig>  mir- <dig> and mir- <dig> with fold change ranging from  <dig>  to  <dig> .
correlation between differentially expressed mirnas and target mrnas
the mrna expression data of matched samples from our previous study  were used for a pairwise correlation analysis  <cit> .
in total,  <dig> mrna probes were differentially expressed  between two breeds.
a pairwise correlation coefficient analysis was computed between  <dig> mirna probes and  <dig> mrna probes.
among the  <dig>  pearson correlation coefficients,  <dig>  negative correlated mirna-mrna pairs were detected at p <  <dig>  and fdr <  <dig>  for correlation between mirna and mrna.
computational target prediction was performed using targetscan and rna hybrid.
after combining the correlation analysis and target prediction results,  <dig> mirna-mrna pairs containing  <dig> genes and  <dig> mature mirna sequences were retained .
the heat map and hierarchical clustering of mirna-mrna pairs based on their expression levels  shown in fig.
<dig> demonstrates an inverse relationship between mrnas and mrna target candidates.
all the target genes were further analyzed with ipa to identify prominent functions and pathways that may contribute to divergent muscle metabolic properties between the two breed types.
target genes assigned to the functional categories related to skeletal and muscular system development and function as well as carbohydrate metabolism were focused on.
the top three canonical pathways for duroc up-regulated target genes were protein ubiquitination pathway, p70s6k signaling and mouse embryonic stem cell pluripotency, while geranylgeranyl diphosphate biosynthesis, phagosome maturation and urate biosynthesis/inosine 5′-phosphate degradation for pinn up-regulated target genes.
a representative mirna-mrna regulatory network of focused biological pathways depicted in fig.
<dig> illustrates a complex relationship and networking of the two biomolecule types.fig.
<dig> heatmaps and hierarchical cluster of mrna and mirna pairs based on least squares  means of differentially expressed mrna and mirna  between breed duroc and pinn
fig.
<dig> differentially expressed mirna-mrna pairs and regulatory network between breed duroc and pinn.
blue indicates mirnas with higher expression in pinn, purple indicates mirnas with higher expression in duroc, red indicates genes with higher expression in duroc and green indicates genes with higher expression in pinn
correlation between mirna expression and phenotypic traits
the expression of  <dig> mirna probes was calculated for the correlation with traits of muscle fiber composition, mitochondrial respiratory activity and metabolic enzyme activity in both duroc and pinn pigs.
in total,  <dig> mirna-phenotype pairs containing  <dig> mirna probes belonged to  <dig> mature mirna sequences were identified at p <  <dig>  shown correlation between mirna and at least one of the  <dig> phenotypes.
the top  <dig> significant mirnas correlated to phenotypes were shown in additional file  <dig> .
table  <dig> showed the top five mirnas significantly correlated to each phenotype.table  <dig> the top five significant mirnas correlated to phenotypes in both duroc and pinn pigs
integration of correlated mirnas, mrnas and phenotypic traits
correlations between gene expression derived from post quality-filtered  <dig>  mrna probes and each phenotypic- trait were calculated for both duroc and pinn pigs.
in total,  <dig>  mrna-phenotype pairs containing  <dig>  mrna probes belonging to  <dig> genes were identified at p <  <dig> .
the top  <dig> significant mrnas correlated with phenotypes are accessible in additional file  <dig> .
pairwise correlation coefficient analysis was then performed between the identified  <dig> mirna probes and  <dig>  mrna probes which correlated with at least one of the  <dig> phenotypes.
after combining with the target prediction results,  <dig>  mirna-mrna pairs containing  <dig> genes and  <dig> mirnas  were identified to correlate with at least one phenotype.
the top ten mirna-mrna pairs for each phenotype were shown in fig.
<dig> and additional file  <dig> .fig.
<dig> regulatory network of mirna-mrna associated with muscle fiber composition, mitochondrial respiratory activity, metabolic enzyme activity and adenine nucleotide concentration for breed duroc and pinn.
genes were colored with orange while mirnas with grey
correlation relationship between mitochondrial and ups related genes
from all identified top  <dig> mirna-mrna pairs associated with each phenotypic trait , the expression correlation between  <dig> selected nuclear-encoded mitochondrial-related genes and  <dig> selected ups-related genes were calculated.
in table  <dig>  mitochondria related genes: atp synthase, mitochondrial f <dig> complex, subunit e , malic enzymes  <dig> , mitochondrial carrier  <dig> , cytochrome p <dig>  family  <dig>  subfamily a, peptide  <dig> , kinesin family member  <dig> binding protein , prohibitin  <dig> , pyruvate dehydrogenase alpha  <dig>  and ubiquinol-cytochrome c reductase complex chaperone  were significantly correlated with at least one of the six ups-related genes: sus scrofa similar to e <dig> ubiquitin-protein ligase hecw <dig> , ubiquitin specific peptidase  <dig> , ubiquitin family domain containing  <dig> , mitochondrial ubiquitin ligase activator of nfkb 1-like , amyloid beta  precursor protein  and heat shock 70 kda protein  <dig>  at p <  <dig> .table  <dig> correlation between mitochondrial and ups related gene expressions
atp5i
me3
mtch1
cyp24a1
kiaa1279
phb2
pdha1
uqcc
hecw2
r
p
<dig>
<dig>
<dig>
<dig>
usp24
r
p
<dig>
<dig>
ubfd1
r
p
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
mul1
r
p
<dig>
<dig>
app
r
p
<dig>
<dig>
hspa2
r
p
<dig>
<dig>
<dig>
<dig>
<dig>
numbers in bold are p-values < significance threshold of  <dig>
qrt-pcr validation
the expression of ssc-mir-24-3p, ssc-mir-30a-5p, ssc-mir- <dig> and ssc-mir- <dig> were random selected for validation by qrt-pcr.
the correlation coefficient between qpcr and michip data ranged from  <dig>   to  <dig>  , suggesting a good concordance between michip and qpcr results, as shown in fig.
<dig> fig.
<dig> qpcr validation of michip results for four micrornas: ssc-mir-24-3p, ssc-mir-30a-5p, ssc-mir- <dig> and ssc-mir- <dig>  plot between qpcr  and michip  for each mirna.
the corresponding correlation coefficient  and p values are shown
discussion
duroc and pietrain are divergent for different muscle characteristics and meat quality.
duroc pigs are fattier and their skeletal muscle contains larger amount sto fibers which are generally associated with higher oxidative enzyme activities, mitochondrial respiratory activity and the storage of lipid to improve the tenderness and juiciness of the meat, whereas pietrain pigs are more muscular and their muscles are more lean and contain higher percentage of ftg fibers which are associated with higher glycolytic enzyme activities .
the higher percentage of ftg fibers may result in lower capillarization, insufficient delivery of oxygen and glycogen depletion, which ultimately lead to dry, firm, and dark meat  <cit> .
since micrornas have been known as critical regulators in energy metabolism of skeletal muscle, in the present study, mirna expression profiling and functional analysis may shed light on mirna-based epigenetic regulatory mechanism of muscle fiber type and metabolic enzyme activity and hence may be translated to improvement of meat quality.
the mirna profile of porcine skeletal muscle has been investigated in several studies and the mirna-mrna networks are constructed .
hou et al.
analyzes the mirna and mrna profile for landrace  and tongcheng  pigs  <cit> .
in both the study of hou et al.
and our present study, biological process of muscle development was shown in the identified differentially expressed mirna-mrna network and it highlights the importance of the regulatory role of mirnas in diverging porcine muscle development between obese-type and lean-type pigs.
further, some mirnas such as mir- <dig>  mir- <dig> and mir- <dig> were identified in the network of both studies.
tang et al.
investigated the mirna and mrna expressions of skeletal muscle for landrace and tongcheng pigs at  <dig>   <dig> and 90 days to explore prenatal muscle development while jing et al.
construct differentially expressed mirna-mrna network between different residual feed intake in pigs  <cit> .
roles of the differentially expressed mirnas and their target genes in divergent muscle characteristics
to understand the regulatory role of mirnas that may contribute to phenotypic variation of the skeletal muscle in different pig breeds, differentially expressed mirnas and mrnas between the two breeds were integrated, so that the key target genes regulated by key mirnas were identified.
functional analysis results showed that the target genes were significantly enriched for various muscle related bio-functions suggesting biological related results rather than random noise.
the top canonical pathway for duroc-up regulated genes was protein ubiquitination pathway containing genes of usp <dig>  usp <dig>  usp <dig>  usp <dig>  usp <dig>  ubr <dig>  dnajc <dig>  fbxw <dig>  birc <dig>  usp <dig>  ube3a, birc <dig>  the atp-dependent ubiquitin is mediated by ubiquitin-activating enzyme e <dig>  specific ubiquitin-conjugating-enzyme e <dig> and ubiquitin protein e <dig> to promote protein degradation via the 26 s proteasome and has implications on meat quality and muscle atrophy .
in our study, usp <dig> and usp <dig> were predicted as a direct target of mir- <dig> and mir- <dig> respectively.
mir-133b could influence both major apoptosis pathways and wound healing  <cit> , and most importantly, the polymorphisms in the porcine mir206/mir133b cluster are proposed as a genetic factor affecting muscle fibers and meat quality traits  <cit> .
mir-311-3p belongs to mir- <dig> mirna family and its loss of function can cause defects in energy metabolism and deregulation of nutritional homeostasis-associated genes  <cit> .
mir- <dig> has been discovered as a negative regulator of adipogenesis in adipose tissue-derived stromal cells by directly targeting the 3′utr of e2f <dig>  <cit> .
this is in line with our study that the expression level of mir- <dig> is higher in duroc than pinn with average fold change more than  <dig>  other duroc up-regulated genes such as cmya <dig>  ar, rb <dig> and bmpr <dig> in functional category of skeletal and muscular system development and function were regulated by mir- <dig>  mir- <dig> and mir- <dig> etc.
cardiomyopathy associated  <dig>  also called trim <dig> belongs to the tripartite motif super family of proteins .
its interaction with both m-band titin and calpain  <dig> suggests its relevance to limb-girdle muscular dystrophies  <cit> .
one snp  of cmya <dig> is significantly associated with water loss and intramuscular fat, which proposes the porcine cmya <dig> as a potential candidate gene for meat quality  <cit> .
androgen receptor  is a steroid-hormone activated transcriptional factor.
the androgen-ar signaling pathway promotes the slow-twitch muscle fiber formation in skeletal muscle by increasing the expression of slow-twitch-specific genes and suppressing the fast-twitch-specific genes  <cit> .. moreover, retinoblastoma  <dig>  has been identified to be related to marbling trait in cattle via gene co-expression analysis  <cit> .
bone morphogenetic protein receptor type ii  encodes a member of the bone morphogenetic protein receptor family of transmembrane serine/threonine kinases.
it is essential for bmp signaling and may be involved in the regulation of adipogenesis and hence in obesity  <cit> .
all the above duroc-up regulated genes with their corresponding down regulated mirnas could contribute to the higher amount of oxidative muscle fibers and fat content.
on the other hand, pinn-up regulated genes smad <dig> and pfkfb <dig> were regulated by mir- <dig>  these two genes are involved in not only muscle fiber specificity but also the promotion of glycolysis in skeletal muscle.
smad family member  <dig>  promotes muscle atrophy in vivo by regulating atrogin- <dig>  mtor and protein synthesis  <cit> .
it suppresses the expression of peroxisome proliferator-activated receptor γ coactivator 1-α   <cit> , which is a master coordinator to control mitochondrial biogenesis and drives the formation of slow-twitch muscle fibers  <cit> .
fructose - <dig>  6-biphosphatase  <dig>  can promote glycolysis by controlling the level of fructose  <dig>  6- bisphosphate, which is an allosteric activator of phosphofructokinase   <cit> .
collectively, we demonstrate that differential mirnas and target gene candidates assemble regulatory networks that may fine-tune the expression of genes within the pathways and shape the related phenotypes among pig breeds.
phenotype-correlated micrornas link to various processes in energy metabolism
numerous studies have revealed the critical roles of mirnas in skeletal muscle development and the irregular mirna expression contributes to various muscular disorders  <cit> .
table  <dig> shows the top five mirnas significantly correlated to each phenotypic trait including muscle fiber type, mitochondrial respiratory activity  and metabolic enzyme activities for duroc and pinn pigs.
both positive and negative correlations were included to discover any potential links between mirnas and phenotypes.
sto, fto, and ftg fibers were three major muscle fiber types in pigs which are strong associated with muscle metabolic activities and meat quality.
mirnas such as mir- <dig>  mir- <dig>  mir- <dig> and mir- <dig> showed highly significant correlation with muscle fiber types.
mir- <dig> and mir- <dig> play a dominant role in the specification of muscle fiber identity by inducing the type i fiber program via targeting sox <dig> which acts as a repressor of slow-twitch genes .
mir-130b belongs to mir- <dig> family directly targets gene pgc-1α, which is a master regulator for mitochondrial biogenesis and its activation promotes the slow, oxidative myogenic program in mice and drives the formation of slow-twitch muscle fibers in cultured muscle cells  <cit> .
in c2c <dig> cells, mir-130b could modulate cellular atp levels by targeting electron transport chain subunits ndufb <dig> and cox6a <dig>  <cit> .
as previously mentioned, mir- <dig> has been discovered as a negative regulator of adipogenesis in adipose tissue-derived stromal cells by directly targeting the 3′utr of e2f <dig>  <cit> .
since lipids are stored mainly in sto muscle fibers to improve tenderness and juiciness of the meat  <cit> , and muscles containing more sto fibers are associated with higher oxidative enzyme activities and mitochondrial respiratory activity  <cit> , it is expected that mir- <dig> shows high correlation with sto muscle fiber type and mitochondrial respiratory activity including state  <dig> pyruvate and state  <dig> succinate.
mir- <dig> family members haven been demonstrated to control calcium signaling by directly inhibiting a ca2+ transporter trpc <dig> etc.
<cit> .
considering the crucial regulatory role of calcium signaling in mitochondrial atp production, it was not unexpected that mir- <dig> was associated with mitochondrial respiratory activity including state <dig> pyruvate and state <dig> succinate and tca involved cs enzyme activity.
furthermore, mir- <dig> directly targets prdm <dig> to promote fast muscle formation since prdm <dig> regulates muscle fiber differentiation by repressing the transcription factor sox <dig> which acts as a repressor of slow-twitch-specific gene expression  <cit> .
all these reinforced the association between sto fiber, mitochondrial respiratory activity and fat content.
mir- <dig> was highly correlated to mitochondrial state  <dig> and state  <dig> respiration rate and complex iv activity.
mir-196a displays a tissue-specific expression pattern in porcine and plays a role in porcine adipose development via inducing the expression of adipocyte specific markers, lipid accumulation and triglyceride content  <cit> .
mir- <dig> was significantly correlated to rci succinate.
mir-542-3p directly targets bone morphogenetic protein  <dig> , which induces differentiation of adipose derived mesenchymal stem cells into brown adipocytes and increases mitochondrial activity in mature brown adipocytes .
mir- <dig> was correlated to complex ii activity in the present study.
mir- <dig> family is abundantly expressed in cardiac and skeletal muscle.
it post-transcriptionally represses the expression of genes in antioxidant network and thus influences susceptibility to cardiac oxidative stress of mir- <dig> transgenic mice  <cit> .
moreover, it was proposed as a candidate gene associated with muscle fiber type composition  <cit> .
mir- <dig>  mir- <dig> and mir- <dig> were identified to correlate with gp activity.
both mir- <dig> and mir- <dig> could directly target and suppress the expression of insulin-like growth factor  <dig> receptor  whereas mir- <dig> regulates insulin synthesis at its mrna level .
since insulin and insulin-like growth factor system are crucial for normal glucose homeostasis  <cit> , it is likely that mir- <dig> and mir- <dig> could play a role in glucose metabolism via igf-1r and insulin.
mir- <dig>  mir- <dig> and mir- <dig> were highly correlated to the concentration of adp and atp in muscle cells.
mir- <dig> and mir- <dig> regulate the expression of oxidative phosphorylation  machinery including complex iv subunits cox <dig>  coxiv and atp synthase subunits atp5g <dig> correspondingly  <cit> .
the fully assemble of oxphos system could directly contribute to the atp production.
furthermore, the modulation of cellular atp levels by mir-15b was supported by other work as well  <cit> .
overall, our results and previous reports functionally link mirnas to muscle fiber specificity, mitochondrial respiration, adipogenesis, glucose metabolism and atp production and further suggest an essential role of mirnas in energy metabolism.
phenotype-correlated mirna-mrna regulatory network link to energy metabolism
based on the identified mirnas that highly correlated with the phenotypes, we further integrated the mirna and mrna expression profiles to identify mirnas regulated genes that influence energy metabolism.
the microrna-mrna regulatory network was constructed using the following criteria: 1) the expressions of both mirnas and target mrnas were correlated to the phenotypical traits 2) the gene expression level was negatively correlated with the expression of its mirna regulator 3) the gene was computationally predicted as a target gene of the corresponding mirna.
it is noteworthy that mir- <dig> together with its target genes bone morphogenetic protein receptor type ii  and insulin receptor substrate  <dig>  were correlated to sto and fto muscle fibers.
mir- <dig> has been documented to be abundant in cardiomyocytes.
it targets the mitochondrial calcium uniporter  and ca2+ transporting atpase  and plays a role in cardiac contractility  <cit> .
in the present study, mir- <dig> was proposed to target both genes bmpr <dig> and irs <dig>  bmpr <dig> encodes a member of the bone morphogenetic protein receptor family of transmembrane serine/threonine kinases.
it is essential for bmp signaling and may be involved in the regulation of adipogenesis and hence in obesity  <cit> .
irs <dig> is a major molecule mediating insulin-signaling pathways.
insulin not only regulates stimulation of protein synthesis and glucose storage  <cit> , but also has effect on mitochondrial function and oxidative capacity of skeletal muscle via increasing the expression level of complex i and complex iv and hence atp production  <cit> .
mir- <dig> and its target gene ubiquitin specific peptidase  <dig>  were correlated to sto fibers, mitochondrial respiratory activity including state  <dig> pyruvate and state  <dig> succinate, and amp concentration in muscle cells, while mir- <dig> and its target gene hecw <dig> were correlated to sto muscle fibers.
usp <dig> belongs to a large family of cysteine proteases that function as deubiquitinating enzymes.
hecw <dig> encodes hect, c <dig> and ww domain containing e <dig> ubiquitin protein ligase  <dig> which is a major component of ubiquitin proteasome system .
ups utilizes atp to promote protein degradation and regulate muscle mass.
accumulated ubiquitin proteins in fast- to slow- transforming rabbit muscle revealed a possible role of ups in muscle fiber specificity  <cit> .
interestingly, mir- <dig> has been discovered as a negative regulator of adipogenesis as described previously  <cit> .
misregulation of mirnas belonging to mir-23a/27a/24- <dig> cluster has been recently associated to hypertrophic cardiomyopathy and skeletal muscle atrophy  <cit> .
mir- <dig> was almost sixfold greater in slow-twitch than in fast-twitch muscle in vivo.
it posttranscriptionally regulates fast-specific myostatin  expression, which mature mrna level is sixfold greater in fast vs slow muscle  <cit> .
in this study, mir- <dig> was identified to be associated with pfk activity – a rate-limiting enzyme in glycolysis and potentially target cytochrome p <dig>  family  <dig>  subfamily a, polypeptide  <dig> , which catalyzes the side-chain oxidation of vitamin d  <cit> .
the vitamin d pathway has the suppressive effect on brown adipocyte differentiation and mitochondrial respiration  <cit> .
mir- <dig> and its predicted targets atp5i, me <dig>  mtch <dig> and cpt <dig> were highly correlated to adp and atp concentration in present study.
mir- <dig> modulates mitochondrial function, decreases cox <dig> expression and activates the generation of reactive oxygen species   <cit> .
atp5i encodes atp synthase mitochondrial f <dig> complex subunit e and it is required for the full assembly of the atp synthase and atp production  <cit> .
me <dig> encodes mitochondrial nadp -dependent malic enzyme  <dig>  the regulation of human mitochondrial nadp -dependent malic enzyme by atp and fumarate may be crucial for the metabolism of glutamine for energy production  <cit> .
mtch <dig> and cpt <dig> encode for mitochondrial carrier  <dig> and carnitine palmitoyltransferase  <dig> respectively.
mtch <dig> also known as presenilin 1-associated protein  which acts as a proapoptotic mitochondrial protein induce apoptosis independent of the proapoptotic proteins bax and bak  <cit> .
the two isoforms of mtch <dig> share two proapoptotic domains and multiple internal signals for import into the mitochondrial outer membrane  <cit> .
fatty acid is a major energy source for the heat and skeletal muscle.
cpt <dig> together with cpt <dig> are involved in beta-oxidation of long chain fatty acids in the mitochondria  <cit> .
altogether suggests that mir- <dig> and its target genes atp5i, me <dig>  mtch <dig> and cpt <dig> are likely to be involved in atp production, apoptosis and beta-oxidation of long fatty acids in mitochondria.
we have demonstrated that correlation relationship between mirna and target mrna can be used to functionally link to phenotypes of interest such as muscle fiber type specification, mitochondrial respiration activity and metabolic enzymes related to atp production.
crosstalk between mitochondria and ups in skeletal muscle
up to now we have shown that the identified mirna-mrna networks are linked to muscle fiber types, oxidative enzyme activities and atp generation.
some of these target genes are involved in mitochondrial function and ups.
an interesting finding was the significant correlation between mitochondrial and ups related gene expressions.
more specifically, mitochondrial related genes including atp5i, me <dig>  mtch <dig>  cyp24a kiaa <dig>  phb <dig>  pdha <dig> and uqcc were highly correlated to at least one of the ups-related genes including hecw <dig>  usp <dig>  ubfd <dig>  mul <dig>  app and hspa <dig>
the tightly interdependent relationship between mitochondria and ups system has been described in many age-related diseases such as alzheimer’s and parkinson’s disease .
our present study revealed a link between these two systems at the level of gene expression under normal state, since all the investigated animals were healthy.
in other words, both mitochondria and ups might contribute to energy metabolism of skeletal muscle via fine-turning the gene expression by mirnas under physiological conditions.
the hect, c <dig> and ww domain containing e <dig> ubiquitin protein ligase  <dig> , ubiquitin specific peptidase  <dig> , ubiquitin family domain containing  <dig>  and mitochondrial ubiquitin ligase activator of nfkb 1-like , amyloid beta precursor protein  and heat shock 70 kda protein  <dig>  are either the major components or associated with the ups system  <cit> .
those genes targeted by several mirnas including mir- <dig>  mir- <dig>  mir- <dig>  mir- <dig>  mir- <dig>  mir- <dig> and mir- <dig> may be the cause of fluctuation of the ups degradation for ubiquitin proteasome-dependent molecules  <cit>  such as transcriptional coactivator pgc-1α which acts as a master regulator for mitochondrial biogenesis, to control mitochondrial gene expression indirectly.
on the other side, mir- <dig> and mir- <dig> targeted genes of atp synthase mitochondrial f <dig> complex subunit e , pyruvate dehydrogenase alpha  <dig>  and ubiquinol-cytochrome c reductase complex chaperone   <cit>  could affect cellular atp generation, followed by influencing the atp-dependent ups system  <cit> .
however, further detail information of the interaction between mitochondria and ups still remains elusive.
CONCLUSIONS
in this study, we modelled the mirna-mrna regulatory networks related to muscle fiber type, metabolic enzyme activity and atp production using the correlation information between expressed mirnas and target mrnas as well as muscular phenotypic measurements of duroc and pinn pigs.
these complex networks may contribute to the muscle phenotypic variations by fine-tuning the expression of genes.
altogether, the results provide an insight into a regulatory role of mirnas in muscular energy metabolisms.
abbreviations
mirna, microrna; sto, slow-twitch-oxidative; fto, fast-twitch-oxidative; ftg, fast-twitch-glycolytic; pinn: pietrain; ups, ubiquitin proteasome system.
additional files
additional file 1: definition and measurement of phenotypical traits.
additional file 2: the primer sequences in qrt-pcr.
additional file 3: differentially expressed mirnas in porcine ld muscle between duroc and pietrain breed types at adult stage.
additional file 4: different expressed mirna-mrna pairs between duroc and pinn after target prediction.
additional file 5: individual mirnas and their correlated phenotypes for duroc and pinn.
additional file 6: individual mrnas and their correlated phenotypes for duroc and pinn.
additional file 7: phenotype correlated mirna-mrna pairs after target prediction for duroc and pinn.

 micrornas are endogenous small non-coding rnas ~22 nt in length that play critical roles in diverse biological processes via epigenetic regulation of gene expression. precursor mirnas are initially generated in nucleus and processed into an approximately 70 nt long stem-loop structure. it is then exported to cytoplasm and processed by dicer to generate mirna/mirna duplexes. one strand of which is incorporated with agonaute to form rna-induced silencing complex that targets mrnas via base-pair complementary, typically to their 3′ untranslated regions or cds and downregulates gene expression by either degradation of mrna or repression of translation, while the other strand is usually discarded. skeletal muscle is highly metabolically active and valuable for meat-producing animals. slow-twitch-oxidative, fast-twitch-oxidative, and fast-twitch-glycolytic fiber were the three major muscle fiber types in pigs. muscle fibers have strong association with muscle metabolic activities and meat quality such as tenderness, juiciness and color. muscle containing a high proportion of oxidative fibers is often associated with higher fat content, oxidative enzyme activities and mitochondrial density; a high ratio of ftg fibers is associated with high glycolytic enzyme activities. previous research has identified several mirnas are associated with meat quality such as mir- mir- and mir- etc. in porcine skeletal muscle. the polymorphisms in the porcine mir- and mir- are proposed as a genetic factor affecting muscle fibers and meat quality traits. since the critical roles of mirnas such as myogenesis, adipogenesis and muscle development have been discovered in pig skeletal muscle, the understanding of the mirna regulation in metabolic properties of skeletal muscle fibers could be key to improvement of meat quality. mir- and mir- could regulate the gene expression of oxidative phosphorylation machinery including complex iv subunits cox coxiv and atp synthase subunits atp5g. mir-15a and mir-15b modulate the cellular atp levels. mir- regulates the fatty acid oxidation and mitochondrial biogenesis through targeting peroxisome proliferator-activated receptor-gamma coactivator-1alpha. with those mirnas identified, the lack of a comprehensive and systematic mirna profiling associated with energy metabolism of skeletal muscle remains unraveled. our previous research on muscle transcriptional profile has revealed numerous biological pathways significantly associated with muscle fiber type, mitochondrial respiratory activity and metabolic enzymes. it is of interesting to further investigate how mirnas are involved in energy metabolism by fine-turning gene expression. in the present study, the mirna transcriptome profiling of longissimus muscle samples obtained 24 h before slaughter of two pig breeds duroc and pietrain exhibiting divergent meat quality and muscle phenotypes provided a comprehensive insight into the discovery of mirnas associated with muscle fiber, mitochondrial respiratory activity and metabolic enzyme activity. muscle of duroc pigs contains higher percentage of sto fibers, mitochondrial respiratory activity and higher fat content to improve the tenderness and juiciness of the meat. in comparison, pinn pigs are more muscular and favorable for meat industry. their skeletal muscles are leaner and contain more ftg fibers. hence duroc and pietrain pigs are great models to study energy metabolism of skeletal muscle. the mirna and mrna expression profile was then integrated based on their pairwise correlations and computational target prediction to construct the regulatory mirna-mrna networks which could potentially affect metabolic properties of skeletal muscle and hence meat quality. the illumination of mirna-based regulatory metabolism could enrich our knowledge of the roles of mirnas in achieving phenotypic diversity of skeletal muscle in different breeds. methods animals and sample collection the experiment and muscle biopsy collection were approved and authorized by the german and european animal welfare regulations for animal husbandry, transport, and slaughter. all experimental procedures, including animal care and tissue sample collection, followed guidelines for safeguarding and good scientific practice in accordance with the german law of animal protection, officially authorized by the animal care committee and authorities. as previously described, duroc and pietrain pigs were raised until 180 days of age. to rule out the effects of the maglinant hyperthermia syndrome locus, only muscle samples from mhs-negative genotype pigs were investigated. muscle biopsies were collected from five female and male pigs of each breed for phenotypic measurements. biopsies were collected 24 h before slaughter from the longissimus muscle between 13th and 14th thoracic vertebrae. muscle samples were frozen in liquid nitrogen and stored at −80 °c until analysis. rna isolation small rnas were isolated and enriched from longissimus muscle biopsies using a mireasy mini kit and an rneasy minelute cleanup kit according to manufacturer’s protocols. the quality and quantity of small rna were assessed with an agilent bioanalyzer using an agilent small rna kit. microrna microarray analysis the affymetrix genechip mirna array was used to determine the mirna expression profile of the lm 24 h ante mortem of duroc and pinn pigs. it is comprised of entries representing hairpin precursor, total probe set for detection of most mirna from species. 200 ng of small rna were used in sample preparation using a flashtag biotin hsr rna labeling kit for affymetrix genechip mirna arrays. each of the labeled rna was then hybridized for 16 h to an affymetrix genechip mirna array according to manufacturer’s recommendation, washed and stained in the affymetrix fluidics station and scanned on the affymetrix g gene array scanner. expression console software was used for robust multichip average normalization and the detection of present mirnas by applying the dabg algorithm. further filtering was done by excluding probes that were present in less than 70 % of the samples within each breed and annotated mirnas that had sequence greater than or equal to 30 nt in length. three thousand five hundred eighty seven probes passed the quality filtering and were used for further analysis. the availability of expression data are in the gene expression omnibus public repository with the geo accession number gse80198: gsm2120718-gsm statistics and bioinformatics analyses differential expression analysis for mirna was performed using the anova procedure in jmp genomics breed was treated as a fixed effect. false discovery rate was used to control an error rate of a multiple-hypothesis testing according to benjamin & hochberg. we used our previous microarray-based mrna expression data to integrate with the differentially expressed mirnas and scan for potential target genes. pearson correlation of mirna and mrna expression levels was calculated. both rnahybrid. and targetscan were used to predict targets of mirnas. rnahybrid is computational software that detects the most energetically favorable hybridization sites of a small rna within a large rna. the mirna probe sets were tested with the following parameters: number of hits per target =  and energy cutoff = −25 kcal/mol and maximal internal or bulge loop size per side =  targetscan was used to predict the target gene candidates based on complementarity of the mirna seed sequence and target binding site on the 5′ utr, 3′ utr and protein coding region of the porcine mrna sequences. hafner et al. and chi et al. have found that argonaute-bound target sites in coding sequences are as numerous as those located in 3′ utr in both hek cells and mouse brain. other research suggests that mirna target sites in 3′utr are more efficient at triggering mrnas degradation while cds and 5′utr located sites can effectively repress translation. xu et al. develop novel computational approach for target prediction with sites located along the entire gene sequences to increase the percentage of true positive targets. hence 5′ utr, cds and 3′utr were included in this study to improve the sensitivity of mirna target identification and avoid a substantial number of missing targets. transcripts that negatively correlated with mirna and predicted as potential targets were further passed to functional analysis. ipa software was used for functional analysis of potential mirna target genes. it categorizes genes based on annotated gene functions and statistically tests for over-representation of functional terms within the gene list using fisher’s exact test. the mirna-mrna regulatory networks were visualized using cytoscape v.. quantitative real time pcr for microrna microarray validation four mirnas related to energy metabolism were validated by qpcr of each individual sample. the same samples were used for both qpcr validation and mirna-chips. the cdna was synthesized from 250 ng isolated mirnas using an ncodetm vilotm mirna cdna synthesis kit according to manufacturer’s protocols. real-time pcr was performed using the express sybr greenertm mirna qrt-pcr kit with premixed rox according to manufacturer’s protocols. all measurements were performed in duplicates. the thermal parameters were 50 °c for 2 min, 95 °c for 2 min, followed by 45 cycles of 95 °c for 15 s and 60 °c for 1 min. the universal qpcr primer was provided in the kit and the mirna-specific forward primers were designed for the mirnas of interest. the designed primer sequence information is accessible in additional file geometric mean of the 5s and u transcription levels was used as an internal standard to normalize the mirna expression value. correlation coefficient analysis between the michip and qpcr data was calculated using sas.  differentially expressed mirnas between duroc and pietrain out of probes quality-filtered probes, probes belonged to unique mature mirna sequences were differentially expressed between duroc and pinn using anova on jmp genomics among these, probes belonged to mature mirna sequences were up-regulated in duroc pigs, while probes belonged to different mature mirna sequences were up-regulated in pinn pigs. among these, mir- mir- and mir- were the top three upregulated mirnas in duroc pigs with fold change ranging from to. whereas the top three upregulated mirnas in pinn were mir- mir- and mir- with fold change ranging from to. correlation between differentially expressed mirnas and target mrnas the mrna expression data of matched samples from our previous study were used for a pairwise correlation analysis. in total, mrna probes were differentially expressed between two breeds. a pairwise correlation coefficient analysis was computed between mirna probes and mrna probes. among the pearson correlation coefficients, negative correlated mirna-mrna pairs were detected at p <  and fdr <  for correlation between mirna and mrna. computational target prediction was performed using targetscan and rna hybrid. after combining the correlation analysis and target prediction results, mirna-mrna pairs containing genes and mature mirna sequences were retained. the heat map and hierarchical clustering of mirna-mrna pairs based on their expression levels shown in fig.  demonstrates an inverse relationship between mrnas and mrna target candidates. all the target genes were further analyzed with ipa to identify prominent functions and pathways that may contribute to divergent muscle metabolic properties between the two breed types. target genes assigned to the functional categories related to skeletal and muscular system development and function as well as carbohydrate metabolism were focused on. the top three canonical pathways for duroc up-regulated target genes were protein ubiquitination pathway, p70s6k signaling and mouse embryonic stem cell pluripotency, while geranylgeranyl diphosphate biosynthesis, phagosome maturation and urate biosynthesis/inosine 5′-phosphate degradation for pinn up-regulated target genes. a representative mirna-mrna regulatory network of focused biological pathways depicted in fig.  illustrates a complex relationship and networking of the two biomolecule types.fig. heatmaps and hierarchical cluster of mrna and mirna pairs based on least squares means of differentially expressed mrna and mirna between breed duroc and pinn fig. differentially expressed mirna-mrna pairs and regulatory network between breed duroc and pinn. blue indicates mirnas with higher expression in pinn, purple indicates mirnas with higher expression in duroc, red indicates genes with higher expression in duroc and green indicates genes with higher expression in pinn correlation between mirna expression and phenotypic traits the expression of mirna probes was calculated for the correlation with traits of muscle fiber composition, mitochondrial respiratory activity and metabolic enzyme activity in both duroc and pinn pigs. in total, mirna-phenotype pairs containing mirna probes belonged to mature mirna sequences were identified at p <  shown correlation between mirna and at least one of the phenotypes. the top significant mirnas correlated to phenotypes were shown in additional file. table  showed the top five mirnas significantly correlated to each phenotype.table the top five significant mirnas correlated to phenotypes in both duroc and pinn pigs integration of correlated mirnas, mrnas and phenotypic traits correlations between gene expression derived from post quality-filtered mrna probes and each phenotypic- trait were calculated for both duroc and pinn pigs. in total, mrna-phenotype pairs containing mrna probes belonging to genes were identified at p < . the top significant mrnas correlated with phenotypes are accessible in additional file. pairwise correlation coefficient analysis was then performed between the identified mirna probes and mrna probes which correlated with at least one of the phenotypes. after combining with the target prediction results, mirna-mrna pairs containing genes and mirnas were identified to correlate with at least one phenotype. the top ten mirna-mrna pairs for each phenotype were shown in fig.  and additional file.fig. regulatory network of mirna-mrna associated with muscle fiber composition, mitochondrial respiratory activity, metabolic enzyme activity and adenine nucleotide concentration for breed duroc and pinn. genes were colored with orange while mirnas with grey correlation relationship between mitochondrial and ups related genes from all identified top mirna-mrna pairs associated with each phenotypic trait, the expression correlation between selected nuclear-encoded mitochondrial-related genes and selected ups-related genes were calculated. in table  mitochondria related genes: atp synthase, mitochondrial f complex, subunit e, malic enzymes, mitochondrial carrier, cytochrome p family subfamily a, peptide, kinesin family member binding protein, prohibitin, pyruvate dehydrogenase alpha and ubiquinol-cytochrome c reductase complex chaperone were significantly correlated with at least one of the six ups-related genes: sus scrofa similar to e ubiquitin-protein ligase hecw, ubiquitin specific peptidase, ubiquitin family domain containing, mitochondrial ubiquitin ligase activator of nfkb 1-like, amyloid beta precursor protein and heat shock 70 kda protein at p < .table correlation between mitochondrial and ups related gene expressions atp5i me3 mtch1 cyp24a1 kiaa1279 phb2 pdha1 uqcc hecw2 r p usp24 r p ubfd1 r p mul1 r p app r p hspa2 r p numbers in bold are p-values < significance threshold of qrt-pcr validation the expression of ssc-mir-24-3p, ssc-mir-30a-5p, ssc-mir- and ssc-mir- were random selected for validation by qrt-pcr. the correlation coefficient between qpcr and michip data ranged from to, suggesting a good concordance between michip and qpcr results, as shown in fig.  fig. qpcr validation of michip results for four micrornas: ssc-mir-24-3p, ssc-mir-30a-5p, ssc-mir- and ssc-mir- plot between qpcr and michip for each mirna. the corresponding correlation coefficient and p values are shown discussion duroc and pietrain are divergent for different muscle characteristics and meat quality. duroc pigs are fattier and their skeletal muscle contains larger amount sto fibers which are generally associated with higher oxidative enzyme activities, mitochondrial respiratory activity and the storage of lipid to improve the tenderness and juiciness of the meat, whereas pietrain pigs are more muscular and their muscles are more lean and contain higher percentage of ftg fibers which are associated with higher glycolytic enzyme activities. the higher percentage of ftg fibers may result in lower capillarization, insufficient delivery of oxygen and glycogen depletion, which ultimately lead to dry, firm, and dark meat. since micrornas have been known as critical regulators in energy metabolism of skeletal muscle, in the present study, mirna expression profiling and functional analysis may shed light on mirna-based epigenetic regulatory mechanism of muscle fiber type and metabolic enzyme activity and hence may be translated to improvement of meat quality. the mirna profile of porcine skeletal muscle has been investigated in several studies and the mirna-mrna networks are constructed. hou et al. analyzes the mirna and mrna profile for landrace and tongcheng pigs. in both the study of hou et al. and our present study, biological process of muscle development was shown in the identified differentially expressed mirna-mrna network and it highlights the importance of the regulatory role of mirnas in diverging porcine muscle development between obese-type and lean-type pigs. further, some mirnas such as mir- mir- and mir- were identified in the network of both studies. tang et al. investigated the mirna and mrna expressions of skeletal muscle for landrace and tongcheng pigs at and 90 days to explore prenatal muscle development while jing et al. construct differentially expressed mirna-mrna network between different residual feed intake in pigs. roles of the differentially expressed mirnas and their target genes in divergent muscle characteristics to understand the regulatory role of mirnas that may contribute to phenotypic variation of the skeletal muscle in different pig breeds, differentially expressed mirnas and mrnas between the two breeds were integrated, so that the key target genes regulated by key mirnas were identified. functional analysis results showed that the target genes were significantly enriched for various muscle related bio-functions suggesting biological related results rather than random noise. the top canonical pathway for duroc-up regulated genes was protein ubiquitination pathway containing genes of usp usp usp usp usp ubr dnajc fbxw birc usp ube3a, birc the atp-dependent ubiquitin is mediated by ubiquitin-activating enzyme e specific ubiquitin-conjugating-enzyme e and ubiquitin protein e to promote protein degradation via the 26 s proteasome and has implications on meat quality and muscle atrophy. in our study, usp and usp were predicted as a direct target of mir- and mir- respectively. mir-133b could influence both major apoptosis pathways and wound healing, and most importantly, the polymorphisms in the porcine mir206/mir133b cluster are proposed as a genetic factor affecting muscle fibers and meat quality traits. mir-311-3p belongs to mir- mirna family and its loss of function can cause defects in energy metabolism and deregulation of nutritional homeostasis-associated genes. mir- has been discovered as a negative regulator of adipogenesis in adipose tissue-derived stromal cells by directly targeting the 3′utr of e2f. this is in line with our study that the expression level of mir- is higher in duroc than pinn with average fold change more than other duroc up-regulated genes such as cmya ar, rb and bmpr in functional category of skeletal and muscular system development and function were regulated by mir- mir- and mir- etc. cardiomyopathy associated also called trim belongs to the tripartite motif super family of proteins. its interaction with both m-band titin and calpain suggests its relevance to limb-girdle muscular dystrophies. one snp of cmya is significantly associated with water loss and intramuscular fat, which proposes the porcine cmya as a potential candidate gene for meat quality. androgen receptor is a steroid-hormone activated transcriptional factor. the androgen-ar signaling pathway promotes the slow-twitch muscle fiber formation in skeletal muscle by increasing the expression of slow-twitch-specific genes and suppressing the fast-twitch-specific genes.. moreover, retinoblastoma has been identified to be related to marbling trait in cattle via gene co-expression analysis. bone morphogenetic protein receptor type ii encodes a member of the bone morphogenetic protein receptor family of transmembrane serine/threonine kinases. it is essential for bmp signaling and may be involved in the regulation of adipogenesis and hence in obesity. all the above duroc-up regulated genes with their corresponding down regulated mirnas could contribute to the higher amount of oxidative muscle fibers and fat content. on the other hand, pinn-up regulated genes smad and pfkfb were regulated by mir- these two genes are involved in not only muscle fiber specificity but also the promotion of glycolysis in skeletal muscle. smad family member promotes muscle atrophy in vivo by regulating atrogin- mtor and protein synthesis. it suppresses the expression of peroxisome proliferator-activated receptor γ coactivator 1-α, which is a master coordinator to control mitochondrial biogenesis and drives the formation of slow-twitch muscle fibers. fructose - 6-biphosphatase can promote glycolysis by controlling the level of fructose 6- bisphosphate, which is an allosteric activator of phosphofructokinase. collectively, we demonstrate that differential mirnas and target gene candidates assemble regulatory networks that may fine-tune the expression of genes within the pathways and shape the related phenotypes among pig breeds. phenotype-correlated micrornas link to various processes in energy metabolism numerous studies have revealed the critical roles of mirnas in skeletal muscle development and the irregular mirna expression contributes to various muscular disorders. table  shows the top five mirnas significantly correlated to each phenotypic trait including muscle fiber type, mitochondrial respiratory activity and metabolic enzyme activities for duroc and pinn pigs. both positive and negative correlations were included to discover any potential links between mirnas and phenotypes. sto, fto, and ftg fibers were three major muscle fiber types in pigs which are strong associated with muscle metabolic activities and meat quality. mirnas such as mir- mir- mir- and mir- showed highly significant correlation with muscle fiber types. mir- and mir- play a dominant role in the specification of muscle fiber identity by inducing the type i fiber program via targeting sox which acts as a repressor of slow-twitch genes. mir-130b belongs to mir- family directly targets gene pgc-1α, which is a master regulator for mitochondrial biogenesis and its activation promotes the slow, oxidative myogenic program in mice and drives the formation of slow-twitch muscle fibers in cultured muscle cells. in c2c cells, mir-130b could modulate cellular atp levels by targeting electron transport chain subunits ndufb and cox6a. as previously mentioned, mir- has been discovered as a negative regulator of adipogenesis in adipose tissue-derived stromal cells by directly targeting the 3′utr of e2f. since lipids are stored mainly in sto muscle fibers to improve tenderness and juiciness of the meat, and muscles containing more sto fibers are associated with higher oxidative enzyme activities and mitochondrial respiratory activity, it is expected that mir- shows high correlation with sto muscle fiber type and mitochondrial respiratory activity including state pyruvate and state succinate. mir- family members haven been demonstrated to control calcium signaling by directly inhibiting a ca2+ transporter trpc etc.. considering the crucial regulatory role of calcium signaling in mitochondrial atp production, it was not unexpected that mir- was associated with mitochondrial respiratory activity including state pyruvate and state succinate and tca involved cs enzyme activity. furthermore, mir- directly targets prdm to promote fast muscle formation since prdm regulates muscle fiber differentiation by repressing the transcription factor sox which acts as a repressor of slow-twitch-specific gene expression. all these reinforced the association between sto fiber, mitochondrial respiratory activity and fat content. mir- was highly correlated to mitochondrial state and state respiration rate and complex iv activity. mir-196a displays a tissue-specific expression pattern in porcine and plays a role in porcine adipose development via inducing the expression of adipocyte specific markers, lipid accumulation and triglyceride content. mir- was significantly correlated to rci succinate. mir-542-3p directly targets bone morphogenetic protein, which induces differentiation of adipose derived mesenchymal stem cells into brown adipocytes and increases mitochondrial activity in mature brown adipocytes. mir- was correlated to complex ii activity in the present study. mir- family is abundantly expressed in cardiac and skeletal muscle. it post-transcriptionally represses the expression of genes in antioxidant network and thus influences susceptibility to cardiac oxidative stress of mir- transgenic mice. moreover, it was proposed as a candidate gene associated with muscle fiber type composition. mir- mir- and mir- were identified to correlate with gp activity. both mir- and mir- could directly target and suppress the expression of insulin-like growth factor receptor whereas mir- regulates insulin synthesis at its mrna level. since insulin and insulin-like growth factor system are crucial for normal glucose homeostasis, it is likely that mir- and mir- could play a role in glucose metabolism via igf-1r and insulin. mir- mir- and mir- were highly correlated to the concentration of adp and atp in muscle cells. mir- and mir- regulate the expression of oxidative phosphorylation machinery including complex iv subunits cox coxiv and atp synthase subunits atp5g correspondingly. the fully assemble of oxphos system could directly contribute to the atp production. furthermore, the modulation of cellular atp levels by mir-15b was supported by other work as well. overall, our results and previous reports functionally link mirnas to muscle fiber specificity, mitochondrial respiration, adipogenesis, glucose metabolism and atp production and further suggest an essential role of mirnas in energy metabolism. phenotype-correlated mirna-mrna regulatory network link to energy metabolism based on the identified mirnas that highly correlated with the phenotypes, we further integrated the mirna and mrna expression profiles to identify mirnas regulated genes that influence energy metabolism. the microrna-mrna regulatory network was constructed using the following criteria: 1) the expressions of both mirnas and target mrnas were correlated to the phenotypical traits 2) the gene expression level was negatively correlated with the expression of its mirna regulator 3) the gene was computationally predicted as a target gene of the corresponding mirna. it is noteworthy that mir- together with its target genes bone morphogenetic protein receptor type ii and insulin receptor substrate were correlated to sto and fto muscle fibers. mir- has been documented to be abundant in cardiomyocytes. it targets the mitochondrial calcium uniporter and ca2+ transporting atpase and plays a role in cardiac contractility. in the present study, mir- was proposed to target both genes bmpr and irs bmpr encodes a member of the bone morphogenetic protein receptor family of transmembrane serine/threonine kinases. it is essential for bmp signaling and may be involved in the regulation of adipogenesis and hence in obesity. irs is a major molecule mediating insulin-signaling pathways. insulin not only regulates stimulation of protein synthesis and glucose storage, but also has effect on mitochondrial function and oxidative capacity of skeletal muscle via increasing the expression level of complex i and complex iv and hence atp production. mir- and its target gene ubiquitin specific peptidase were correlated to sto fibers, mitochondrial respiratory activity including state pyruvate and state succinate, and amp concentration in muscle cells, while mir- and its target gene hecw were correlated to sto muscle fibers. usp belongs to a large family of cysteine proteases that function as deubiquitinating enzymes. hecw encodes hect, c and ww domain containing e ubiquitin protein ligase which is a major component of ubiquitin proteasome system. ups utilizes atp to promote protein degradation and regulate muscle mass. accumulated ubiquitin proteins in fast- to slow- transforming rabbit muscle revealed a possible role of ups in muscle fiber specificity. interestingly, mir- has been discovered as a negative regulator of adipogenesis as described previously. misregulation of mirnas belonging to mir-23a/27a/24- cluster has been recently associated to hypertrophic cardiomyopathy and skeletal muscle atrophy. mir- was almost sixfold greater in slow-twitch than in fast-twitch muscle in vivo. it posttranscriptionally regulates fast-specific myostatin expression, which mature mrna level is sixfold greater in fast vs slow muscle. in this study, mir- was identified to be associated with pfk activity – a rate-limiting enzyme in glycolysis and potentially target cytochrome p family subfamily a, polypeptide, which catalyzes the side-chain oxidation of vitamin d. the vitamin d pathway has the suppressive effect on brown adipocyte differentiation and mitochondrial respiration. mir- and its predicted targets atp5i, me mtch and cpt were highly correlated to adp and atp concentration in present study. mir- modulates mitochondrial function, decreases cox expression and activates the generation of reactive oxygen species. atp5i encodes atp synthase mitochondrial f complex subunit e and it is required for the full assembly of the atp synthase and atp production. me encodes mitochondrial nadp -dependent malic enzyme the regulation of human mitochondrial nadp -dependent malic enzyme by atp and fumarate may be crucial for the metabolism of glutamine for energy production. mtch and cpt encode for mitochondrial carrier and carnitine palmitoyltransferase respectively. mtch also known as presenilin 1-associated protein which acts as a proapoptotic mitochondrial protein induce apoptosis independent of the proapoptotic proteins bax and bak. the two isoforms of mtch share two proapoptotic domains and multiple internal signals for import into the mitochondrial outer membrane. fatty acid is a major energy source for the heat and skeletal muscle. cpt together with cpt are involved in beta-oxidation of long chain fatty acids in the mitochondria. altogether suggests that mir- and its target genes atp5i, me mtch and cpt are likely to be involved in atp production, apoptosis and beta-oxidation of long fatty acids in mitochondria. we have demonstrated that correlation relationship between mirna and target mrna can be used to functionally link to phenotypes of interest such as muscle fiber type specification, mitochondrial respiration activity and metabolic enzymes related to atp production. crosstalk between mitochondria and ups in skeletal muscle up to now we have shown that the identified mirna-mrna networks are linked to muscle fiber types, oxidative enzyme activities and atp generation. some of these target genes are involved in mitochondrial function and ups. an interesting finding was the significant correlation between mitochondrial and ups related gene expressions. more specifically, mitochondrial related genes including atp5i, me mtch cyp24a kiaa phb pdha and uqcc were highly correlated to at least one of the ups-related genes including hecw usp ubfd mul app and hspa the tightly interdependent relationship between mitochondria and ups system has been described in many age-related diseases such as alzheimer’s and parkinson’s disease. our present study revealed a link between these two systems at the level of gene expression under normal state, since all the investigated animals were healthy. in other words, both mitochondria and ups might contribute to energy metabolism of skeletal muscle via fine-turning the gene expression by mirnas under physiological conditions. the hect, c and ww domain containing e ubiquitin protein ligase, ubiquitin specific peptidase, ubiquitin family domain containing and mitochondrial ubiquitin ligase activator of nfkb 1-like, amyloid beta precursor protein and heat shock 70 kda protein are either the major components or associated with the ups system. those genes targeted by several mirnas including mir- mir- mir- mir- mir- mir- and mir- may be the cause of fluctuation of the ups degradation for ubiquitin proteasome-dependent molecules such as transcriptional coactivator pgc-1α which acts as a master regulator for mitochondrial biogenesis, to control mitochondrial gene expression indirectly. on the other side, mir- and mir- targeted genes of atp synthase mitochondrial f complex subunit e, pyruvate dehydrogenase alpha and ubiquinol-cytochrome c reductase complex chaperone could affect cellular atp generation, followed by influencing the atp-dependent ups system. however, further detail information of the interaction between mitochondria and ups still remains elusive.  in this study, we modelled the mirna-mrna regulatory networks related to muscle fiber type, metabolic enzyme activity and atp production using the correlation information between expressed mirnas and target mrnas as well as muscular phenotypic measurements of duroc and pinn pigs. these complex networks may contribute to the muscle phenotypic variations by fine-tuning the expression of genes. altogether, the results provide an insight into a regulatory role of mirnas in muscular energy metabolisms. abbreviations mirna, microrna; sto, slow-twitch-oxidative; fto, fast-twitch-oxidative; ftg, fast-twitch-glycolytic; pinn: pietrain; ups, ubiquitin proteasome system. additional files additional file 1: definition and measurement of phenotypical traits. additional file 2: the primer sequences in qrt-pcr. additional file 3: differentially expressed mirnas in porcine ld muscle between duroc and pietrain breed types at adult stage. additional file 4: different expressed mirna-mrna pairs between duroc and pinn after target prediction. additional file 5: individual mirnas and their correlated phenotypes for duroc and pinn. additional file 6: individual mrnas and their correlated phenotypes for duroc and pinn. additional file 7: phenotype correlated mirna-mrna pairs after target prediction for duroc and pinn. 
in the current study, mirna expression profile was investigated in longissimus muscle biopsies of malignant hyperthermia syndrome-negative duroc and pietrain pigs with distinct muscle metabolic properties in order to explore the regulatory role of mirnas related to mitochondrial respiratory activity and metabolic enzyme activity in skeletal muscle. a comparative analysis of the mirna expression profile between duroc and pietrain pigs was performed, followed by integration with mrna profiles based on their pairwise correlation and computational target prediction. next, we analyzed the correlation between individual mirnas and phenotypical traits including muscle fiber type, mitochondrial respiratory activity, metabolic enzyme activity and adenosine phosphate concentrations, and constructed the regulatory mirna-mrna networks associated with energy metabolism. these complex networks relate to muscle fiber type, metabolic enzyme activity and atp production and may contribute to divergent muscle phenotypes by fine-tuning the expression of genes. musclemitochondrial respiratory activitymirna-mrna networkmeasurements of mitochondrial respiratory activity and metabolic enzyme activity were supported by the german research foundation issue-copyright-statement© the author 2016 altogether, the results provide an insight into a regulatory role of mirnas in muscular energy metabolisms and may have an implication on meat quality and production. in the identified mirna-mrna network, a tight relationship between mitochondrial and ubiquitin proteasome system at the level of gene expression was observed. the identified target genes were enriched in protein ubiquitination pathway, stem cell pluripotency and geranylgeranyl diphosphate biosynthesis, as well as skeletal and muscular system development. we assembled mirna-mrna regulatory networks based on divergent muscle properties between different pig breeds and further with the correlation analysis of expressed genes and phenotypic measurements. micrornas are small non-coding rnas that play critical roles in diverse biological processes via regulation of gene expression including in skeletal muscles. it revealed a link between these two systems contributing to energy metabolism of skeletal muscle under physiological conditions. it is noteworthy that mir- targeting bmpr and irs mir- targeting usp mir- targeting hecw and mir- targeting atp5i, me mtch and cpt were highly associated with slow-twitch oxidative fibers, fast-twitch oxidative fibers, adp and atp concentration suggesting an essential role of the mirna-mrna regulatory networking in modulating the mitochondrial energy expenditure in the porcine muscle.
 micrornas are small non-coding rnas that play critical roles in diverse biological processes via regulation of gene expression including in skeletal muscles. in the current study, mirna expression profile was investigated in longissimus muscle biopsies of malignant hyperthermia syndrome-negative duroc and pietrain pigs with distinct muscle metabolic properties in order to explore the regulatory role of mirnas related to mitochondrial respiratory activity and metabolic enzyme activity in skeletal muscle.  a comparative analysis of the mirna expression profile between duroc and pietrain pigs was performed, followed by integration with mrna profiles based on their pairwise correlation and computational target prediction. the identified target genes were enriched in protein ubiquitination pathway, stem cell pluripotency and geranylgeranyl diphosphate biosynthesis, as well as skeletal and muscular system development. next, we analyzed the correlation between individual mirnas and phenotypical traits including muscle fiber type, mitochondrial respiratory activity, metabolic enzyme activity and adenosine phosphate concentrations, and constructed the regulatory mirna-mrna networks associated with energy metabolism. it is noteworthy that mir- targeting bmpr and irs mir- targeting usp mir- targeting hecw and mir- targeting atp5i, me mtch and cpt were highly associated with slow-twitch oxidative fibers, fast-twitch oxidative fibers, adp and atp concentration suggesting an essential role of the mirna-mrna regulatory networking in modulating the mitochondrial energy expenditure in the porcine muscle. in the identified mirna-mrna network, a tight relationship between mitochondrial and ubiquitin proteasome system at the level of gene expression was observed. it revealed a link between these two systems contributing to energy metabolism of skeletal muscle under physiological conditions.  we assembled mirna-mrna regulatory networks based on divergent muscle properties between different pig breeds and further with the correlation analysis of expressed genes and phenotypic measurements. these complex networks relate to muscle fiber type, metabolic enzyme activity and atp production and may contribute to divergent muscle phenotypes by fine-tuning the expression of genes. altogether, the results provide an insight into a regulatory role of mirnas in muscular energy metabolisms and may have an implication on meat quality and production. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords musclemitochondrial respiratory activitymirna-mrna networkmeasurements of mitochondrial respiratory activity and metabolic enzyme activity were supported by the german research foundation issue-copyright-statement© the author 2016 