Deciphering the microRNA transcriptome of skeletal muscle during porcine development

Animal ethics statement

All animals used in this study were farmed according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and approved by the Institutional Animal Care and Use Committee in the College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China under permit No. DKY-B20121406.

Preparation of experimental animals and tissues

Fifteen female Jinhua pigs, which were grouped into five stages (three replicates for each stage): prenatal embryonic day 90 (E90), postnatal 0 d (natal), 30 d (weaning), 180 d (young adult), and 7 y (middle-aged), were humanely sacrificed. The longissimus dorsi muscle tissues from the same anatomical location were collected from pigs representing each of the five stages. Subsequently, the samples were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction. Total RNA was isolated from each sample using Trizol (Takara) according to the manufacturer’s protocol. The concentration and purity of RNA samples were determined by measuring the A260/A280 ratio using a NanoDrop ND1000 spectrophotometer (Thermo Scientific, Hudson, NH, USA).

Myofiber histology

After sacrifice, all muscle tissues were fixed in 10% neutral buffered formalin solution, embedded in paraffin using a TP1020 semi-enclosed tissue processor (Leica, Wetzlar, Germany), sliced at a thickness of 6 µm using RM2135 rotary microtome (Leica, Wetzlar, Germany), and stained with hematoxylin and eosin (H&E). The mean diameter of muscle cells was calculated as the geometric average of the maximum and minimum diameter, and 100 cells were measured for each sample in randomly selected fields using a TE2000 fluorescence microscope (Nikon, Tokyo, Japan) and Image Pro-Plus 7.0 software.

Small RNA library construction and high-throughput sequencing

We pooled the total RNA for replicates in the E90, 0-d, and 30-d stages. For 180-d and 7-y stages, total RNA of each replicate was individually used for library construction. Therefore, a total of nine libraries were constructed and sequenced using single-end sequencing in 36 nt reads by an Illumina Genome Analyzer II. The bioinformatics pipeline for miRNA discovery and profiling was carried out as previously described, with some improvements (Gambardella et al., 2010). All reads were counted and the identical reads were combined into a single kind. After eliminating adaptor sequences, low-quality tags, sequences smaller than 16 bp, and reads with no insertion, all of the clean tags were annotated and classified by comparison with the non-coding RNAs (rRNA, tRNA, scRNA, snRNA, and snoRNA) in the GenBank (http://www.ncbi.nlm.nih.gov/) and Rfam 9.1 (http://rfam.sanger.ac.uk/) databases. Known miRNAs were identified by comparing our clean tags to mature miRNAs in miRBase 19.0. The remaining non-annotated sRNA sequences were aligned against the porcine genome, and genomic sequences containing the sRNA with 80–100 nt of flanking sequences were used to predict hairpin structures using Mfold (http://mfold.rna.albany.edu). Only the sequences that could be folded into typical hairpin structures and were located in intergenic regions or introns were considered to be miRNA precursor loci of potential novel miRNAs in the porcine genome.

miRNA differential expression and clustering analyses

IDEG6 (Romualdi et al., 2003) was employed to detect DE miRNAs between two libraries. A unique miRNA was considered DE when P < 0.001 was yielded by three statistical tests (Greller and Tobin, R of Stekel and Falciani, and general chi-squared tests) with Bonferroni correction. Hierarchical clustering analysis (HCL) of DE miRNAs was performed with MultiExperiment Viewer (MeV) (Howe et al., 2010).

Prediction and functional annotation of miRNA target genes

The potential targets of a certain miRNA were predicted by PicTar (Krek et al., 2005) and TargetScan Human 6.2 (Lewis, Burge & Bartel, 2005), and the pairwise overlaps of the results from both programs composed the final predicted targets. The predictions were based on human mRNA–miRNA interactions, because of the absence of porcine miRNAs in the current version of abovementioned algorithm. The enriched GO biological process (GO-BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms of predicted target genes were determined using the DAVID bioinformatics resource (Huang, Sherman & Lempicki, 2008).

STEM analysis

STEM analysis was used to visualize expression patterns for DE miRNAs. Each miRNA was assigned to the model profile that its time series most closely matched based on the correlation coefficient. The number of miRNAs assigned to each model profile was then computed. The number of miRNAs expected to be assigned to a profile was estimated by randomly permuting the original time point values, renormalizing the miRNA expression values, assigning miRNAs to their most closely matching model profiles, and repeating this process for a large number of permutations. The average number of miRNAs assigned to a model profile over all permutations was used as the estimate of the expected number of miRNAs assigned to the profile. The statistical significance of the number of miRNAs assigned to each profile compared with the expected number was also computed (Ernst & Bar-Joseph, 2006).

qRT-PCR validation

The expression changes of nine selected miRNAs were validated using an EvaGreen-based High-Specificity miRNA qRT-PCR Detection Kit (Stratagene, La Jolla, CA, USA) on the CFX96™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The q-PCR validations were carried out on three biological replicates. Information on the primer pairs used is available in Table S1. Three endogenous control genes (U6 snRNA, 18S rRNA, and 5S rRNA) were used in this assay. The ΔΔCt method was used to determine the expression level differences between surveyed samples.

Phenotypic measurements

In this study, we investigated muscle fiber diameter and body weight change during the postnatal muscle development stage. As shown in Figs. S1A–S1C, we found that muscle fiber diameter and body weight were significantly increased (P < 0.001, one-way analysis of variance) from 0 d to 7 y, which is consistent with results from previous studies in pigs (Dwyer, Fletcher & Stickland, 1993; Fiedler et al., 1997). Notably, we found that muscle fiber diameter and body weight growth rate significantly slowed with age (P < 0.05) (Figs. 1A and 1B). This age-related decrease of muscle fiber growth rate is likely due to decreases in protein synthesis and increases in protein degradation (Jackman & Kandarian, 2004). These phenotypic differences indicated possible underlying intrinsic differences in molecular mechanisms.

Summary of deep-sequencing data

To investigate miRNA expression changes during pig muscle development, nine small RNA libraries were constructed using deep sequencing. In total, 85.84 million raw reads (with redundancy) were obtained. More than 56.67% reads (ranging from 56.67% to 81.25%) in each library passed the Adapter (ADT) dimmer & length, junk, mRNA, RFam, and Repbase filters and were considered high-quality reads (Fig. S2). Of these high-quality reads (Fig. S3), the majority (90.35 ± 1.34%, n = 9) of the small RNAs were approximately 21–24 nt in length. More than half of the reads were 22 nt in length (57.03 ± 7.53%, n = 9), followed by 23 nt (16.40 ± 4.34%, n = 9), 21 nt (12.37 ± 5.26%, n = 9), and 24 nt (4.55 ± 1.49%, n = 9), which are typical lengths of Dicer-processed products (Bartel, 2009; Berezikov, 2011). This result indicates that small RNA sequencing is a reliable approach to generate miRNA reads for further analysis.

miRNA transcriptome profiles during skeletal muscle development

The high-quality reads (sequences) were compared against the pig genome and divided into three groups (Table 1): porcine known miRNAs, porcine conserved miRNAs, and porcine putative new miRNAs. Of the porcine known miRNAs, 404 corresponded to 233 porcine known pre-miRNAs in miRBase 19.0 (http://www.mirbase.org/) and could be mapped to the pig genome. Specifically, 286 miRNAs were known, but 118 were not previously identified porcine miRNAs (Table S2). Of the porcine conserved miRNAs, 101 corresponded to 71 other known mammalian species except for porcine pre-miRNAs in miRBase 19.0, which were further mapped to the pig genome (Table S3). Of the porcine putative new miRNAs, 431 corresponded to 325 candidate pre-miRNAs that were not mapped to mammalian pre-miRNAs in miRBase 19.0, but predicted RNA hairpins derived from the pig genome, and the extended genome sequences form hairpins, which were labeled PN (Table S4). Therefore, 629 pre-miRNAs that encode 936 mature miRNAs were identified. Because there are distinct pre-miRNAs and genomic loci that code identical mature sequences, 903 unique miRNAs were finally identified (Table S5).

To uncover the possible roles of these abundant miRNAs in each stage, we ranked the miRNAs by expression level. The profile of miRNA expression abundance was different across different development stages. As shown in Fig. S4, the most abundant miRNAs (top 10 unique miRNAs) account for more than 50.58% counts of all 903 unique miRNAs. Particularly, the counts of top 10 miRNAs occupied about 90% of the total counts in the 180-d and 7-y libraries, which is significantly higher than that of the E90, 0-d, and 30-d libraries. Furthermore, the set of the top 10 unique miRNAs over the five muscle development stages corresponded to 16 kinds of unique miRNAs (Fig. 2). Of these miRNAs, six (miR-133a-1/-2-3p, let-7a-1/-2-5p, miR-27b-3p, miR-26a-5p, miR-1-3p, and let-7f-1/-2-5p) were shared by all five stages and were closely related to myogenesis, cell growth, myocyte proliferation, and cell apoptosis. For example, miR-133a is a muscle-enriched miRNA that inhibits proliferation of progenitor cells and promotes myogenesis by targeting SRF (McCarthy & Esser, 2007). miR-26a promotes myogenesis in C2C12 cells via post-transcriptionally repressing Ezh2, which is a known suppressor of skeletal muscle cell differentiation (Wong & Tellam, 2008). These two miRNAs are most abundant during the fast-growing stage of pig muscle. Moreover, miR-27a/b, a potential regulator of myogenesis, could induce skeletal muscle hypertrophy by down-regulating myostatin, an inhibitor of myogenesis (Huang et al., 2012; Sharma et al., 2014) and miR-27b inhibition leads to more proliferation and delays the onset of differentiation (Crist et al., 2009). miR-1, a muscle-specific microRNA, promotes cell apoptosis by targeting Bcl-2 (Tang et al., 2009), and could target heat shock protein 70 (HSP70) in the development of muscle atrophy (Kukreti et al., 2013). Therefore, miR-27b promotes myogenesis and proliferation, whereas miR-1 inhibits these processes and induces apoptosis. miR-27b showed increasing expression levels during development, and miR-1 showed the opposite trend. These results may indicate fluctuation of miRNA regulation during muscle development.

In addition, 10 other miRNAs (miR-378-1/-2-3p, miR-127-3p, miR-191-5p, miR-486-2-5p, miR-143-3p, miR-10a-5p, miR-148a-3p, miR-99a-5p, miR-30e-5p, and miR-199a-1/-2-5p) (Fig. 2) in the set of the top 10 most highly expressed unique miRNAs over the five muscle development stages are related to various cell proliferation, myogenesis, and apoptosis responses. For example, miR-378 promotes myogenesis in pigs through regulation of bone morphogenetic protein 2 (BMP2) and mitogen-activated protein kinase 1 (MAPK1), especially in fiber formation in both the fetal and newborn periods (Hou et al., 2012). In our study, miR-378 showed the highest expression level in the E90 stage. Though Hou et al. found that the expression level of miR-378 increased at 65 and 90 dpc and peaked at postnatal day 0, these results both suggested that miRNA-378 was a new candidate miRNA for myogenesis in pigs (Hou et al., 2012; Qin et al., 2013). These results supports the idea that it has an important role in fetal muscle development. Another miRNA, miR-148a, belong to the top ten expressed miRNAs only in E90 period, which was in line with previous finding that the average abundance of this miRNA before birth was eight times higher than that in postnatal (Qin et al., 2013). miR-486 is directly controlled by SRF, myocardin-related transcription factor-A (MRTF-A), and MyoD, and regulates muscle growth and homeostasis (Small et al., 2010). Moreover, miR-191 could promote cell proliferation and suppress apoptosis in MGC803 cells (Shi et al., 2011). These miRNAs were within the top 10 expressed miRNAs in 180-d and 7-y pig muscle, but not in relatively early stages, which indicates that they have regulatory roles in muscle function maintenance after sexual maturity. In contrast, miR-199a is a potential regulator of myogenesis through suppression of WNT-signaling factors and regulation of muscle growth and homeostasis (Alexander et al., 2013); its expression was higher in the E90, 0-d, and 30-d stages than in the 180-d and 7-y stages, which indicates that miR-199a is an important myogenic regulatory factor in early fetal and newborn muscle development.

Previous study performed by Qin et al. indicated that most of the highly expressed miRNAs in porcine skeletal muscle such as miR-1 and miR-133 will be more functional. Besides the miR-1 and miR133a that we found highly expressed in all the stages (McCarthy & Esser, 2007), miR-26a also showed abundant expression (Huang, Sherman & Lempicki, 2008; Huang et al., 2008). Taken together, these results indicate that highly abundant miRNAs across various development stages might play important roles in maintenance of vital muscular physiological functions and regulation of skeletal muscle development. Nonetheless, these highly expressed miRNAs also showed dynamic performance across different stages, which indicates typical spatial and temporal-specific expression of miRNAs during muscle development. Undoubtedly, dysregulation of such miRNAs would perturb muscular homeostasis and may even initiate pathogenesis (Eisenberg et al., 2007; McCarthy & Esser, 2007).

Differentially expressed (DE) miRNA during muscle development stages

To determine if miRNA expression significantly differed between libraries, IDEG6 (Romualdi et al., 2003) was employed to normalize calculations between the high-quality sequences in the nine libraries. A unique miRNA was considered differentially expressed when a general chi-squared test produced P-values <0.001 (Li et al., 2011). By applying this criterion, we identified 505/903 miRNAs (55.92%) that were differentially expressed during the muscle development stages (Table S6). A previous study demonstrated that miRNAs function in a dose-dependent manner (Carlsbecker et al., 2010); thus, only the relatively more abundant miRNAs (>1,000 read counts) were included in the subsequent analysis. Therefore, 167 of 505 (33.07%) differentially expressed (DE) miRNAs were retained (>1,000 read counts).

To further characterize variability in miRNA expression profiling, we performed a hierarchical clustering analysis (HCL) based on the DE miRNAs. As shown in Fig. 3, the miRNA expression profiles showed two major clusters: one that included E90, 0-d, and 30-d stages, and one that included 180-d and 7-y stages. This clustering pattern may correspond to the different growth rate of the muscle fibers in different stages (fast-growth stages: E90, 0-d, and 30-d; slow-growth stages: 180-d and 7-y) (Dwyer, Fletcher & Stickland, 1993; Rehfeldt et al., 2000). This significant variation in miRNA expression may be responsible, at least in part, for the phenotypic differences (Fig. 1A). Meanwhile, three biological replicates at 180-d and 7-y were highly correlated (180-d: average r = 0.9532; 7-y: average r = 0.9709), which indicates high experimental reliability and good reproducibility of replicates (Fig. 3).

Furthermore, nine representative DE miRNAs (miR-133a-5p, miR-181a-1-3p, miR-499-5p, miR-320-3p, miR-24-1-3p, miR-214-3p, let-7g-5p, miR-23a-3p, and miR-10b-3p) were chosen for validation by the stem–loop real time (RT)-PCR-based method using three independent samples. The results were consistent with our sequencing result (Pearson r = 0.875 ± 0.119, n = 9; Fig. 4), which indicates that deep sequencing does allow for the successful discovery of miRNAs from pig muscle with high accuracy and efficiency.

Distinct miRNA expression patterns during muscle development

To visually illustrate the expression pattern of the miRNAs during different development stages, a short time-series expression miner (STEM) analysis was performed for the DE miRNAs. As shown in Fig. 5A, four model profiles had a larger than expected number of genes assigned to the cluster and were significant (P < 0.05). For expression pattern 1, expression level decreased across all five stages. For expression pattern 2, expression level increased from E90 to 0 d, then deceased from 0 d to 180 d, and remained stable from 180 d to 7 y. For expression pattern 3, expression level decreased from E90 to 180 d, then deceased from 180 d to 7 y. For expression pattern 4, expression level kept increasing throughout the five stages. The miRNA list and normalized value of expression by STEM software for each pattern was shown in Table S7. These four profiles may indicate the different roles of corresponding miRNAs; therefore, we performed target prediction and Gene Ontology (GO) analysis. As shown in Fig. 5B, the target genes of highly expressed miRNAs in E90 and 0 d (patterns 1, 2, and 3) were mainly involved in processes related to skeletal muscle development; for example, “cellular component morphogenesis,” “embryonic morphogenesis,” and “cell morphogenesis involved in differentiation.” In mammals, muscle mass is mainly determined by the number and size of muscle fibers (Rehfeldt et al., 2000). The number of muscle fibers is prenatally determined during primary and secondary muscle fiber formation (about 30–60 d post-gestation) (Zhao et al., 2003), whereas postnatal muscle fiber growth is mainly associated with accumulation of muscle-specific proteins (hypertrophy). Therefore, these miRNAs with high expression levels in the embryonic period, which are closely related to cellular morphogenesis and embryonic morphogenesis, may promote differentiation of the muscle fibers in the embryonic period.

In contrast, the target genes of highly expressed miRNAs in the 180-d and 7-y stages (Fig. 5B, pattern 4) were mainly involved in processes such as “apoptosis,” “regulation of apoptosis,” “regulation of cell death,” and “positive regulation of macromolecule biosynthetic process,” which are consistent with the results of previous studies (Dirks & Leeuwenburgh, 2005). Postnatal growth of skeletal muscle mainly occurs through increases in muscle fiber length and girth and is inhibited in relatively older animals (Rehfeldt et al., 2000). Moreover, it has been demonstrated that muscle atrophy induced in middle-aged or older animals is associated with higher amounts of cell apoptosis and death (Dirks & Leeuwenburgh, 2005). This may be the primary reason these miRNAs that play vital roles in apoptosis are highly expressed in middle-aged stage. Particularly, miR-1 showed significantly increasing expression level in 180d and 7y stages (Fig. 5A and Table S7). By using both C2C12 myotubes and dex-induced muscular atrophy mouse models, Kukreti et al. (2013) indicated that miR-1 is a muscle-specific microRNA and has a role in promoting muscle atrophy. On the other hand, McCarthy & Esser (2007) also revealed that miR-1 decreased during mouse skeletal muscle hypertrophy. Our result was consistent with the findings in mouse. Hence, miR-1 could be a potential modulator in regulating porcine postnatal skeletal muscle development.

These results together indicate that there are obviously disparate roles of miRNAs in different development stages of skeletal muscle, such as enhanced levels of miRNAs related to differentiation and morphogenesis in the embryonic stage. Conversely, there were increasing levels of apoptosis-related miRNAs relatively later in muscle development, especially in the middle-aged stage.