Integrative Analysis of Porcine microRNAome during Skeletal Muscle Development

Pig is an important agricultural animal for meat production and provides a valuable model for many human diseases. Functional studies have demonstrated that microRNAs (miRNAs) play critical roles in almost all aspects of skeletal muscle development and disease pathogenesis. To investigate the miRNAs involved in regulating different periods of skeletal muscle development, we herein performed a comprehensive research for porcine microRNAome (miRNAome) during 10 skeletal muscle developmental stages including 35, 49, 63, 77, 91 dpc (days post coitum) and 2, 28, 90, 120, 180 dpn (days postnatal) using Solexa sequencing technology. Our results extend the repertoire of pig miRNAome to 247 known miRNAs processed from 210 pre-miRNAs and 297 candidate novel miRNAs through comparison with known miRNAs in the miRBase. Expression analysis of the 15 most abundant miRNAs in every library indicated that functional miRNAome may be smaller and tend to be highly expressed. A series of muscle-related miRNAs summarized in our study present different patterns between myofibers formation phase and muscle maturation phase, providing valuable reference for investigation of functional miRNAs during skeletal muscle development. Analysis of temporal profiles of miRNA expression identifies 18 novel candidate myogenic miRNAs in pig, which might provide new insight into regulation mechanism mediated by miRNAs underlying muscle development.


Introduction
Pork is an important source of protein for human consumption, particularly in China where it tops the list of meat consumed. In addition, the pig can be considered as a suitable model for human diseases because it is similar with human in size, physiology, pathology and genomics [1]. Muscle development is a complex process, including not only embryonic myogenesis but adult myofiber maturation. Most research suggested that porcine primary myofibers mainly formed at 35 to 64 dpc and secondary myofibers at 54 to 90 dpc [2]. However, Swatland et al. pointed out that fetal myofiber hyperplasia ceased at slightly later at 70 days gestation and the increase in number of myofibers seen in a transverse section was caused by growth in length of myofibers [3]. Similarly, Zhao et al. discovered that myogenesis was almost completed before 77 dpc [4], suggesting that the muscle fiber formed at 35 to 77 dpc. Other research described that there existed a third generation of fibers around birth in pigs and might be part of mechanisms leading to the larger muscle mass [5,6]. The total fiber number (TFN) does not increase after birth and is fixed at the number formed during embryonic life. Myofibers were classified according to their contractile and metabolic properties as type 1, type 2 (2A and 2B) and an intermediate type 2X. From late gestation through the first four postnatal weeks, myofibers underwent a process of maturation from type 1 to type 2 [4,6]. Hypertrophy is another process of skeletal muscle maturation, involving an increase in muscle mass and in the size, associated with increased myofibrillar protein content. Signaling pathways such as insulin-like growth factor I (IGF-I) and Ca 2+ /calmodulin-dependent transcriptional pathways have been demonstrated to govern skeletal muscle hypertrophy and atrophy [7]. Understanding the complex mechanism underlying porcine muscle development might contribute to swine industry as well as human muscular growth and diseases.
Increasing evidence suggests that miRNAs serve as biological regulators by mediating gene expression. It was increasingly clear that miRNAs involve almost all aspects of skeletal muscle development, including cell migration, proliferation, differentiation and apoptosis [8,9]. Despite the significant role of miRNAs in the regulation of myogenesis, the miRBase (release 18.0) lists only 257 distinct mature miRNA sequences in pigs, significantly fewer than that in human. It's worth noting that many miRNAs are expressed in a tissuespecific or stage-specific manner [10], and the bestcharacterized muscle-specific miRNAs (myomiRs [11]) are miR-1, miR-206 and miR-133 families which specifically expressed in cardiac and skeletal muscles. Dicer was known to be essential for miRNA biogenesis, whose elimination in the skeletal muscles could lead to prenatal death, muscle mass reduction or myofiber malformation [12], an outcome highlights the significance of miRNAs in muscle development. Currently, scientists tend to study the molecular mechanism of swine skeletal muscle development through analysis of miRNA transcriptome profiles [13,14], discovering and identifying more miRNAs based on bioinformatics analysis. However, it is still not enough for investigation in mapping a nearly complete porcine miRNAome during skeletal muscle development.
The Landrace, one of the most popular pig breeds in most countries in Europe, North America and Asia, was employed in this investigation. Its miRNAome during skeletal muscle development at ten developmental stages including five prenatal stages (35,49,63,77, 91 dpc, simplified as LR 1, 2, 3, 4, 5) and five postnatal stages (2,28,90,120, 180 dpn, simplified as LR 6,7,8,9,10) using Solexa sequencing. The time points cover almost all morphological and physiological changes of porcine skeletal muscle growth and development [4]. Thus, this study will provide a thorough investigation of miRNAome in porcine skeletal muscle to facilitate a better understanding of their involvement in myogenesis.

Ethics statement
All animal procedures were performed according to guidelines developed by the China Council on Animal Care and protocols were approved by the Animal Care and Use Committee of Guangdong Province, China. The approval ID or permit numbers are SCXK (Guangdong) 2004-0011 and SYXK (Guangdong) 2007-0081.

Sample preparation
Fifteen Landrace (LR) purebred sows with the same genetic background were artificially inseminated with semen from the same purebred boars. The pregnant sows were slaughtered at five prenatal stages (35,49,63,77,91 dpc, days post coitum) while the piglets and adult pigs were slaughtered at five postnatal stages (2, 28, 90, 120, 180 dpn, days postnatal). The longissimus dorsi muscle tissues collected per time point were used as the experimental samples for sequencing. These samples were snap-frozen in liquid nitrogen and stored at -80°C.

Small RNA library construction and Solexa sequencing
Total RNA was extracted using miRNeasy Mini Kit (Cat#217004, QIAGEN, GmBH, Germany) according to the manufacturer's protocol. For each developmental stage, equal quantities of total RNA isolated from three individual pigs were pooled. Total RNA integrity was measured on an Agilent 2100 Bioanalyzer system (Agilent) for quality control.16-35 nt RNA fragments were excised, purified from a PAGE gel, and ligated with 5′ and 3′ adaptors using T4 RNA ligase. Reverse transcription followed by PCR was used to create cDNA constructs based on the small RNA ligated with 3′ and 5′ adapters. Subsequently, the amplified cDNA constructs were purified from agarose gel, in preparation for sequencing analysis using the Illumina Genome Analyzer (Illumina, CA, USA) according to the manufacturer's instructions.

Data analysis
The raw data were processed using Illumina 1G Genome Analyzer Pipeline software and then submitted to data filtration. Clean reads were obtained after filtering low-quality reads and trimming the adaptor sequences. All of the clean reads were initially searched against miRBase (version 18; http:// www.mirbase.org/) to identify known porcine miRNAs. Unmappable reads subsequently were annotated and classified by reference to non-coding RNAs in the Ensemble (ftp:// ftp.ensembl.org/pub/release-69/fasta/sus_scrofa/ncrna/), piRNA (http://pirnabank.ibab.ac.in/) and Rfam (version 10; http://rfam.sanger.ac.uk/) databases. The mappable sequences were achieved and used for further analysis. Meantime, many unannotated sequences that cannot match any above databases were analyzed by miRDeep [15]

Target Prediction and KEGG Orthology Analysis
Based on the sequences of the miRNAs, the targets of DE miRNAs were predicted using the online database miRecords [17] (http://mirecords.biolead.org/) and chosen when they were predicted in at least databases integrated in miRecords. The predicted target genes were subsequently submitted to KOBAS [18] for KEGG Orthology analysis (http://kobas.cbi.pku.edu.cn/ home.do) using KEGG database. Given the P values in KOBAS, pathways with statistically significant values (P<0.05) were chosen. In addition, SPSS was used to divide the pathways into categories by using the fuzzy theories to Kmeans algorithm.

Stem-loop RT-PCR
To validate the sequencing results, nine miRNAs with different expression levels were selected and Stem-loop RT-PCR was performed using the Lightcycler480 (Roche) with SYBR-Green detection (SYBR PrimeScript RT-PCR Kit, TaKaRa Biotechnology Co. Ltd.) according to the manufacturer's instructions. The reverse transcriptase reaction was performed using Reverse Transcription System Kit (Promega) with the stem-loop primers listed in Table S7. Each real-time PCR system contained 5 µl 2×SYBR Premix DimerEraser TM (TaKaRa), 0.3 µl forward and reverse primers respectively, 1 µl template cDNA, and dH 2 O up to the final volume of 10 µl. The reactions were incubated at 95°C for 30 s, followed by 40 cycles of 95°C for 10 s, 60°C for 20 s and 72°C for 10 s. Mir-17-5p was chosen as internal control and all reactions were run in triplicate. The experimental data were analyzed using the 2 -∆∆CT method.

Overview of sequencing data
After trimming of adaptor sequences and removal of low quality reads, 181,224,007 total clean reads were obtained from ten libraries, of which almost 96% matched small RNA databases mentioned in the Materials and Methods, namely annotated reads (Table S1). Saturation plots of ten libraries were created to assess the efficiency and accuracy of deep sequencing for miRNA detection ( Figure S1). Taking the LR1 library as an example, the deeper were sequenced, the more unique small RNAs were found. However, the most of unique small RNAs were unknown, for the number of annotated small RNAs accounted for less than 20%. Furthermore, the number of new sequences observed for known small RNAs and porcine miRNAs (miRBase found) reached a plateau when the number of sequenced reads was 9,500,000, suggesting that the library capacity approached saturation. Similar plots of the other nine libraries were shown in Figure S1. Therefore, the results demonstrated that the deep sequencing data were able to accurately represent the miRNA transcriptome profiles of porcine skeletal muscle.
Of the annotated sequences, the most abundant size class in the small RNA sequences distribution was 22 nt, followed by 21 and 23nt ( Figure 1A and B), consistent with the known 21-23 nt range for mature miRNAs. In addition, the annotated sequences were analyzed referencing the data from miRBase (release 18.0, containing 228 precursors and 257 mature miRNAs), resulting 247 (96%) mature miRNAs and 210 (92%) precursor were identified. Next, the rest annotated sequences were subjected to piwi-interacting RNA (piRNA) database, Rfam (version 10) and non-coding RNA in Ensemble and the results were summarized in Table S1. Most annotated small RNAs were porcine miRNAs (ssc-miRNAs), accounting for 76.7% of the total sequence reads but a less proportion (8.7%) of the unique sequence reads in ten small RNA libraries ( Figure 1C). The results indicated that the majority of small RNAs were annotated miRNAs while other classes added diversity. Together, all of the above results provide confidence that the deep sequencing data were highly enriched for known miRNA sequences, suggesting that the data are reliable for analyzing miRNA expression profiles as well as for predicting novel miRNAs.

Identification of potential novel miRNAs
Deep sequencing is a robust approach for discovering novel miRNAs that are expressed at low levels. To find more potential miRNAs, unannotated sequences longer than 18 nt were searched against porcine genome and analyzed by miDeep. According to criteria for novel miRNAs identification, 297 candidate novel miRNAs were predicted and named ssc-miR-new-N (N=1 ~ 297). These novel miRNAs ranged from 20 to 23 nt in length, among which 22 nt sequences accounted for the most (40%) followed by23 nt (25%), 21 nt (22%) and 20nt (13%), consistent with typical length distribution of mature miRNAs. Total reads of these novel miRNAs counted up from ten libraries were at a relatively low level: 267 putative novel miRNAs were sequenced below 100 times, 27 were below 1000 times (of which 3 were above 500 times), 2 were above 1000 times and only 1 were above 10000 times (Table S2). Hairpin structures of the putative novel miRNAs were then analyzed using RNAfold.

Characterization of miRNAome in different stages of skeletal muscle
As presented in Table S3, 247 identified mature porcine miRNAs have been substituted as 257 mature porcine miRNAs (If the normalized expression of a given miRNA is zero, its expression value will be modified to 0.01) that expressed at different levels, ranging from single reads for the least abundant to millions of reads for the most abundant. Previous study has demonstrated that over 60% of miRNAs detected by deep sequencing had no discernible activity, all of which expressed below 100-1,000 reads per million (RPM) [19]. Thus, the distribution of numbers for normalized miRNAs was summarized in Table S3, from which we could find that the number of miRNAs with mean expression value below 1000 is 197 (123 and 74, respectively), almost accounts for 77% of the total number of miRNAs (48% and 29%, respectively) in the ten libraries. However, the number of miRNAs expressed on average over 1000 RPM is 60 (23%), among which 15 miRNAs (6%) expressed over 10000 RPM. Moreover, we counted the top most abundant 15 miRNAs in each library, accounting for 85% of the total counts of all unique miRNAs on average ( Figure S2), from which we obtained a common 27 miRNAs microRNAome during Skeletal Muscle Development PLOS ONE | www.plosone.org with the highest expression level in ten libraries (Table 1, details see Table S4), suggesting that the majority of abundant miRNAs are from fewer miRNAs [20]. Previous studies have announced that there were more miRNAs with reported myogenic function in muscle besides well-known myomiRs ( Table 2). Among the 27 miRNAs, we found 11 (40%) miRNAs have a relationship with muscle development validated by experiments. Almost all known myomiRs were identified with the most abundance. The most abundant miRNA was ssc-miR-1, which presented by more than 2,100,000 RPM in ten libraries. The predominance of miR-1 is consistent with its well established function during skeletal muscle development [21] and reported role during porcine myogenesis [22]. Similarly, two other myomiRs, miR-133 [21] and miR-206 [23], were highly expressed and ranked the 4 th and 6 th respectively, while two other miRNAs (miR-378 [24,25] and miR-143 [25]) ranked the 2 nd and 3 rd have been identified to participate in the proliferation and differentiation of muscle cells.
To identify whether different biological processes responding to different stages of muscle development, KEGG Orthology analyses were performed using predicted targets of DE miRNAs at 35 to 77 dpc (LR 1-4), 77 dpc to 28 dpn (LR4-7) and 28 to 180 dpn (LR7-10). Figure 3B illustrated that most pathways with reported function in regulating myogenesis were significantly enriched (P<0.01) in the above three developmental stages, including Wnt, MAPK, Neurotrophin, TGF-beta, ErbB, Insulin, Calcium, GnRH, T cell receptor signaling pathways and Focal adhesion as well as Adherens junction. In particular, two well-known signaling pathways that Wnt and MAPK were both highly enriched in all of the three stages. Comparing these pathways among the three stages, Wnt, the most differently expressed signaling pathway, was significantly enriched throughout the muscle development, with the highest expression at 28 to 180 dpn, the second at 77dpc to 28 dpn, and the least at 35 to 77 dpc. However, MAPK signaling pathway expressed at the same level in all stages. Other pathways such as TGF-beta and Calcium signaling  pathways showed more differences in different stages, and Insulin, GnRH and ErbB signaling pathways were just the opposite. Notably, almost all pathways showed s similar expression pattern in two waves of myofiber maturation: myofiber transformation and hypertrophy, suggesting their closer relationship in terms of myogenic regulation. Furthermore, we refined a more detailed KEGG Orthology analysis using predicted targets of DE miRNAs in Table 3 between two time points contributing to describing the different developmental phases of myogenesis. The significant pathways (P<0.05) related to muscle development were listed in Table 4, most of which have been identified according to the most abundant miRNAs ( Figure 3A). Further analysis of these pathways suggested three clusters of enrichment by using SPSS, setting values for statistic significance ( * P<0.05 and ** P<0.01 were assigned 1 and 1.5, respectively, using 0 where no significance was detected) : A) with universal enrichment, including Adherens junction, Focal adhesion, Adipocytokine, GnRH, Insulin, Notch, p53, Phosphatidylinositol, TGF-beta, VEGF and Wnt signaling pathways; B) with more enrichment during prenatal myogenesis, including Axon guidance, Gap junction, Melanogenesis and Neurotrophin signaling pathways; C) with more enrichment during postnatal myogenesis, including Calcium, Chemokine, ErbB, Hedgehog, mTOR, T cell receptor and MAPK signaling pathways (Table 4).

Clustering miRNA expression
We used the Short Time-series Expression Miner (v 1.1, STEM) [51] to cluster temporal profiles of miRNA expression during muscle fiber formation at 35 to 77 doc (LR1-4) as well as fiber further maturation at 77 dpc to 180 dpn (LR4-10). After that a total of 57 out of 87 DE miRNAs (65.5% of 87 DE miRNAs; Table S5) observed at 35 to 77 dpc were significantly clustered into three kinds of expression patterns (11.5% of 26 clusters; Figure 4A), while 112 out of 166 DE miRNAs (67.5% of 166 DE miRNAs; Table S5) observed at 77 dpc to 180 dpn were significantly clustered into two kinds of expression patterns (4% of 50 clusters; Figure 4B, the default P-value was 1E-5). Figure 4C showed miRNA expression profiles for the five clusters found to be significant out of 76 possible clusters successively. The cardinality of each cluster was ranged from 7 DE miRNAs in cluster 3 to 80 in cluster 4 (Table S6). Visual examination of these clusters suggested that down-regulated DE miRNAs were distributed both Cluster 1 and 3 while Cluster 2 included DE miRNAs that showed more gradual increase, which peaked at 63 to 77 dpc during myofiber formation. STEM Similar expression pattern could be associated with functional correlation to a certain extent. Notably, in addition to porcine myomiRs (miR-1, -206, -133 a-3p/a-5p/b), three other miRNAs (miR-128, -208b and -378) reported to be related to muscle development in other mammals were also included in up-regulation Cluster 2, suggesting their roles in the regulation of porcine embryonic myogenesis at 35 to 77 dpc. However, Cluster 5 illustrated that ssc-miR-206 and other four musclerelated miRNAs (ssc-miR-126, -148a/b and -15b) continued to decline while ssc-miR-133b and eleven other muscle-related miRNAs (ssc-miR-125b, -128,-181a/b, -199a, -214, -23a, -24, -424, -503 and -7) in Cluster 4 presented a down and then up trend from 77 dpc to 180 dpn, suggesting their different roles played in adult fiber maturation. Moreover, parts of these miRNA expression profiles during 10 developmental stages were illustrated in Figure 5, providing the possibility for further research into miRNAs that had similar expression profiles with muscle-related miRNAs.

Validation of the Sequencing data
Stem-loop quantitative RT-PCR was applied to validate the sequencing data. Choosing an appropriate set of endogenous control (EC) miRNA genes is crucial, for EC genes are widely used to normalize the miRNA q-PCR data, which are expected to express constantly at all stages of development of one or more tissues. Although U6 snRNA is one of the most widely used internal control, Gu et al. demonstrated that the U6 gene was the least stable gene compared with other candidate EC genes in all tissues (including muscle tissues) comparing with miR-17, -103, -107 and -23a [52]. Hence, the optimal EC gene from above candidates was studied by calculating the standard deviation (Stdev.) of four genes in all samples when the fifth was supposed to be an EC gene ( Figure S3). As a consequence, miR-17-5p was actually the most stable gene (Stdev. = 0.1), while the U6 gene was the least stable gene (Stdev. = 4.9), consistent with the previous study [52]. Stemloop quantitative RT-PCR was then performed on 9 random miRNAs with different expression levels (miR-1, -206, -133a-3p, -133a-5p, -133b, -378, -214, -744 and let-7f) to validate the sequencing data ( Figure 6). Stem-loop RT-PCR primers were presented in Table S7. The Pearson correlation coefficient of the Real-time PCR and the Solexa sequencing was calculated and the r values ranged from 0.84 to 0.95, indicating that there was a high consistency between the two methods.

Discussion
Since pigs are important farm animals and are suitable models for studying human diseases, we herein presented the systematical porcine miRNAome at ten important developmental stages of skeletal muscle using Solexa sequencing. A total of 247 known and 297 potential novel  [14], this study aimed at studying the mechanism of skeletal muscle development in swine through analysis of miRNA transcriptome profiles at designed different developmental stages. In addition, we both focused on the changes in abundance of myomiRs (miR-1, -206 and -133) during swine skeletal muscle development. Contrarily, different developmental stages were used to study muscle development in the two investigations. Moreover, the proliferating satellite cells were employed in the previous work, which might be complement each other and contribute to revealing the expression profiles and biological function of miRNA during swine skeletal muscle development more comprehensively. Unlike the previous work, the present study paid more attention to the most abundant miRNAs, but not limited to the well-studied myomiRs, which have been proven to play important roles during swine skeletal muscle development, consistent with the conclusion that functional miRNAs tend to be highly expressed [19]. The sequencing data was then verified by stem-loop quantitative RT-PCR. Moreover, the abundance of the top 15 miRNAs were found to be the majority of total reads in each library ( Figure S2). KEGG Orthology analysis of the most abundant 25 common miRNAs (excluding let-7 family) from 10 libraries indicated that their involvement in the regulation of canonical myogenesis-related pathways and diseases ( Figure 5). Intriguingly, 11 of the 25 common miRNAs were reported to function as muscle-related miRNAs including the well-known myomiRs (miR-1, -206 and -133 family). Therefore, it is reasonable to hypothesize that the functional miRNA, in terms of mediating target suppression to regulate myogenesis and muscle diseases, tended to express abundantly.
Currently, differently expressed (DE) miRNAs were thought to be closely related to almost all aspects of muscle development. In this study, a total of 183 DE miRNAs between different libraries were identified. Ssc-miR-1 was the most abundant miRNA in ten libraries, consistent with the wellestablished function of miR-1 during skeletal muscle development [21]. Interestingly, the second abundant DE miRNA was miR-378, a new candidate miRNA for myogenesis in pigs by down-regulating porcine BMP2 or MAPK1 [53]. Another DE miRNA (miR-148a), whose average abundance before birth was eight times higher than that in postnatal, might be a part of mechanism implicated in differences between embryonic myogenesis and adult myofiber maturation. Exceptions were some previously reported muscle-related miRNAs such as miR-181c/d-5p, miR-29b/c, miR-221/222 and miR-208, all of which expressed at a very low level (average RPM <100). In addition to highly expressed myomiRs, the expression of muscle-related miRNAs was varied in pigs, suggesting that highly expressed muscle-related miRNAs that have been validated in other species should be the first priority in studying porcine skeletal muscle development, followed by the muscle-related miRNAs expressed at lower levels.  Next, target prediction and KEGG Orthology analysis of stage-specific DE miRNAs might suggest that myofiber formation and maturation were controlled by multiple signal pathways in different ways, suggesting that different molecular regulation mechanism underlined different biological processes. In particular, the two waves of myofiber maturation: myofiber transformation and hypertrophy showed a closer relationship in myogenic regulation, providing a possibility to undergo subsequent STEM clustering by combining them. Furthermore, a more detailed KEGG Orthology analysis of DE miRNAs between two time points identified a list of pathways participating in the regulation of myogenesis and related diseases, suggesting that DE miRNAs were tightly related with biological processes at different stages of muscle development. For example, Wnt signaling pathway involved in embryonic myogenesis and in regulating the homeostasis of adult muscle [45], resulting its universal enrichment from embryonic to adult myofiber maturation. Calcium signaling, though with less attention, was considered potential mediators of postnatal muscle development and hypertrophy [54], giving some indication for more enrichment during postnatal myogenesis. An exception was Notch signaling pathway, though it was classified to the category with universal enrichment, only expressed highly during late embryogenesis. However, previous studies have showed its crucial role not only for embryogenesis but for muscle maintenance and repair in the adult [55]. One reason might explain the difference is that there is still plenty of unknown information about miRNA, resulting insufficiency of its functional annotation. In general, integrated analysis of sequencing results and biological pathway information might add insight into the vital roles played by miRNA during skeletal muscle development. However, comprehensive information about miRNA is still lacking, such as insufficiency of functional analysis, which requires more research focusing on combination of bioinformatic analysis of miRNAome and subsequently experimental validation.
We hypothesized that miRNAs with similar expression patterns might be function-related. STEM clustering results suggested that ssc-miR-378 functioned as a new candidate miRNA for porcine myogenesis because of its expression profile similar to ssc-miR-1 and -133a-3p ( Figure 5A). Interestingly, recently studies have validated that ssc-miR-378 regulated myogenesis by directly targeting the BMP2 or MAPK1 in pig, suggesting that the STEM clustering was reliable for analyzing miRNA expression profiles as well as for predicting of candidate myogenic miRNAs. We focused on either highly expressed DE miRNAs or muscle-related miRNAs to ensure the accuracy of prediction. Consequently, 18 candidate miRNAs were selected, including ssc-miR-378, -127, -128, -411, 23b, -27b, -10a, -140*, -9-1/-2, -148a/b, -126, 542-3p, 30a-5p/d/e-5p and miR-103 ( Figure 5). MiR-127, which is located within a CpG island with little research on myogenesis, showed similar expression pattern to ssc-miR-206, up regulated at 35 to 77 dpc and down regulated at 77 dpc to 180 dpn ( Figure 5B), providing the possibility for further investigation concerning the function of ssc-miR-127 during muscle development. miR-128 was reported to participate in the regulation of adipogenesis, osteogenesis and myogenesis [36] and herein, together with ssc-miR-411, up regulated at 35 to 77 dpc and fluctuated at 77 dpc to 180 dpn ( Figure 5C), might behave in a similar manner as ssc-miR-133b during porcine muscle development. MiR-23b was found to interact with TGF bata signaling by down-regulating Smads in fetal hepatocytes [56], while miR-27b was involved in myogenic differentiation [31] as well as fast-specific and glucocorticoid-dependent myostatin expression [30], both of which were highly expressed at 35 dpc and down-regulated thereafter ( Figure 5D), suggesting their roles played in porcine embryonic myogenesis. Like ssc-miR-23b and 27b, it is reasonable to hypothesize that ssc-miR-10a, -140* and -9-1/-2 might act on earlier embryonic myogenesis, for they highly expressed at 35 dpc then decline dramatically. MiR-148a has been identified as a novel myogenic miRNA that mediated myogenic differentiation via targeting ROCK1 [27], while miR-126 attenuated insulin signaling [57] and governed vascular integrity and angiogenesis [58], suggesting their interactions with signaling pathways were required for muscle normal development and maintenance. Similarly, ssc-miR-148b, -542-3p and -30 family (a-5p/d/e-5p) showed similar expression patterns with ssc-miR-148a and -126, highly expressed and down-regulated at 77 dpc to 180 dpn ( Figure  5E), making it possible that they belong to the candidate myogenic miRNAs. MiR-103/107 family was validated to attenuate miRNA biosynthesis by targeting Dicer [59], regulate insulin sensitivity through down-regulating caveolin-1 [60] and affect cellular migration by modulating CDK5R1 expression [61], allowing us to hypothesize that these miRNA-mediated mechanism may affect myogenesis when the microRNA family, especially ssc-miR-103, down-regulated during both myofiber formation (at 35 to 77 dpc) and maturation (at 77 dpc to 180 dpn) ( Figure 5F). In summary, STEM clustering of miRNA temporal expression contributed to determining candidate myogenic miRNAs in an effective way.
In this study, 297 novel miRNAs were identified. Analysis of the evolutionary conservation of these novel miRNAs revealed that none is conserved in mammals. Therefore, it is reasonable to hypothesize that pig-specific miRNAs may exist. However, it is the first time that more than 200 novel miRNAs were detected from miRNAome of porcine skeletal muscle attributing to high-throughput deep sequencing technology. Although few of these novel miRNAs might play any roles in establishing and maintaining phenotype of muscle tissue during individual development for their extremely low expression, there are still several highly expressed novel miRNAs that need further validation. In particular, a functional validation of novel and species-specific miRNAs is a challenge for understanding the critical roles played by miRNAs during muscle development and for providing the valuable information for pig meat quality improvement. Figure S1. Saturation plots of ten libraries.

Supporting Information
(TIF) Figure S2. Counts characteristics of the unique miRNAs in each library. Starting from the miRNA with the highest counts (x-axis), the black bar represents the accumulative proportion of miRNAs in total counts of each library. The red horizontal  Figure S3. Determination of the optimal endogenous control (EC) pig genes for normalization. Based on standard deviation, the stability of 5 candidate EC genes was measured.
The pairwise difference between all tested genes was compared using t-test, and the significance was labeled with asterisk ( ** p value < 0.001). (TIF) Table S1. Annotations of sequenced small RNAs.