Simultaneous miRNA and mRNA transcriptome profiling of human myoblasts reveals a novel set of myogenic differentiation-associated miRNAs and their target genes

Identification of new myogenesis-related MIRNA (MRMIRNA) in human primary myoblasts

To identify new myogenesis-related miRNAs in human cells, CD56+ myogenic precursors were isolated from skeletal muscle biopsies of six healthy individuals (Additional file 1). These were then induced to differentiate in vitro using serum starvation and the progression of myogenic differentiation was monitored by following the expression of Myogenin (MYOG) and other myogenic differentiation markers including ID1, TNNT1, TNNT2, MYL4, COL15A1, VIM and CAV1 (Figure 1A). The disappearance of proliferating cells and concomitant emergence of myotubes was monitored using Ki67 staining and microscopic observation (Figure 1B and C).

We observed a statistically significant difference of the expression levels of myogenesis-related markers in proliferating or differentiated myogenic cells isolated from different healthy subjects. This difference was considered as a natural heterogeneity of gene expression levels in different subjects that might be explained by their age, sex, physical activity or diet preferences, therefore, all samples analyzed were included in the study.

The expression of 365 different microRNAs was then profiled in human primary myoblasts and myotubes using TaqMan Low Density Array (TLDA) (Additional file 2). In total, 60 miRNAs were found to be differentially expressed during myogenic differentiation with p-values ranging from 3.4x10^-5 to 5x10^-2 (Figure 2). These will be referred to as myogenesis-related miRNAs (MR-miRs). Sequences and corresponding TaqMan probe codes for these microRNAs can be found in the Additional file 3.

The majority (43 out of 60) of MR-miRs were found to be upregulated of which 14 had not been previously reported as involved in the regulation of myogenesis either in humans or in other species (Table 1). Of the 17 MR-miRs that were found to be downregulated during myogenesis, 6 had not been previously associated with myogenic differentiation.

MRNA transcriptome profiling and functional classification of differentially expressed genes

Until now, mRNA expression profiling of human myoblasts has been performed exclusively on primary cells [36]. However, primary skeletal myoblasts have not been enriched for myogenic cells prior to analysis. To avoid analyzing mixed populations of myoblasts and non-myogenic cells, we have used RNA extracted from CD56+ myoblasts purified from muscular tissues. Then, using Agilent transcriptome 44 K microarrays, we identified 6,616 differentially expressed genes of which 3,445 and 3,245 were down- or up-regulated, respectively (Additional file 4). In most cases, a functional class could be assigned to these transcripts using the DAVID Gene Ontology database (Additional file 5) [37, 38].

Predication of myogenesis-related MIRNA target genes

A common practice to boost target gene predictions accuracy is to use several prediction algorithms and then to combine their predictions. However, this approach has been proved inefficient [39]. Therefore, we have used a single algorithm, e.g. RNA22 [40] to predict miRNA targets. RNA22 has been chosen for its low rate of false-positive predictions as compared to other algorithms [39]. For each miRNA, the number of target predictions varied between 261 and 3,058.

According to previous reports, the proportion of correct predictions can be as low as 10% [41, 42] (for recent review see [43]. To increase the accuracy of bioinformatic predictions, we used a simultaneous mRNA/miRNA expression profiling approach described in [42, 44, 45]. Based on the fact that MiRNA regulate gene expression by inhibiting translation or inducing deadenylation of mRNAs followed by their degradation (for review see [46, 47]), we reasoned that the expression levels of mRNA and miRNA should be inversely correlated if one regulates the other. Indeed, for upregulated MR-miRNA, the majority of their predicted target genes was found to be underexpressed. Similarly, most downregulated MR-miRNAs had upregulated predicted targets. Indeed, for every microRNA tested, the distribution of Pearson correlation coefficients of its predicted target genes was shifted towards negative values. This was the case neither for randomly selected genes, nor when the sample correspondence was permuted (Figure 3A and the Additional file 6) which was used to establish a significance threshold of Pearson correlation coefficients. Target genes with correlation coefficients below the threshold had less than 5% chances to be random. We call these genes “supported targets” to underscore the fact that the prediction was supported by mRNA transcriptome profiling (Additional file 7). On average, 11% of microRNA target genes predicted by the RNA22 algorithm were “supported” by mRNA transcriptome profiling (Figure 3B).

We then compared the accuracy of predictions made by RNA22 and TargetScan, another popular miRNA target prediction algorithm. As an indication of accuracy, we have used a percentage of predicted targets that were supported by transcriptome data. For 11 of 17 target genes, RNA22 demonstrated higher percentage of predictions supported by transcriptome data, for 2 more microRNAs the difference in the accuracy of predictions was insignificant (Figure 3C and Additional file 8).

QRT-PCR validation of supported MIRNA target genes

To demonstrate the relevance of our miRNA target predictions, we altered the expression of several microRNAs and measured the expression of their predicted target genes, supported by the transcriptome data. To illustrate the impact of miRNA repression on the expression of their target genes, we have chosen to inhibit the expression of four well-characterized myogenic microRNAs miR-1, -133a, -133b and -206 in human immortalized myoblasts using LNA (Locked Nucleic Acids) antisense microRNAs (anti-miRs). Human immortalized myoblasts were differentiated in vitro and then transfected separately with anti-miRs targeting miR-1, miR-133a, miR-133b and miR-206. miR-1 and miR-206 have the same seed sequence and, therefore, target very similar lists of genes. The same is true for miR-133a and miR-133b. We have then randomly selected 12 and 10 genes predicted to be targeted by miR-1/206 and miR-133a/b respectively and supported by our transcriptome data and tested their expression using qRT-PCR in the cells transfected with corresponding anti-miRs. Repressing miR-1/206 resulted in the upregulation of six of its targets: AP3D1 (adaptor-related proteins complex 3, delta 1 subunit), COL3A1 (collagen, type III, alpha 1), HDAC1 (histone deacetylase 1), PDCD4 (programmed cell death 4), PRKAB2 (protein kinase, AMP-activated beta 2 non-catalytic subunit) and ZNF365 (zinc finger protein 365). These target genes were considered as novel “qRT-PCR validated” target genes of miR-1 and -206 (Figures 4A). Six other target genes supported by our transcriptome data including ATP2B1, C11orf45, FOXP, G2E3, MBOAT7 and TUBB could not be validated via qRT-PCR experimentally, as we did not observe significant changes of their expression following anti-miR transfection (data not shown).

In the case of miR-133a and -133b, immortalized myotube transfection anti-miRs resulted in the upregulation of six of ten randomly selected target genes, supported by our transcriptome data. These novel qRT-PCR validated targets of miR-133a and -133b include ASAP2 (ArfGAP with SH3 domain, ankyrin repeat and PH domain 2), ARFGAP1 (ADP-ribosylation factor GTPase activating protein 1), ATP1A2 (ATPase, Na+/K+ transporting, alpha 2 polypeptide), C7orf51 (chromosome 7 open reading frame 51), RBM15B (RNA binding motif protein 15B) and TMCC2 (transmembrane and coiled-coil domain family 2) (Figures 4B). Four other target genes supported by our transcriptome data including LMCD1, NAT9, PRPF38A and SYNGR2 could not be validated via qRT-PCR, as we did not observed significant changes of their expression following anti-miR transfection (data not shown).

To illustrate the impact of miRNA overexpression on the expression of their target genes, we have selected two microRNAs, miR-30 and -128, that were previously shown to be upregulated during myogenic differentiation [26, 27], but their function in myogenesis is not well-established. Two plasmids coding for the corresponding pre-miRNAs or scrambled sequence were transfected separately into human rhabdomyosarcoma cells, and the expression of their target genes was then tested via qRT-PCR. Ectopic overexpression of miR-128 led to the downregulation of three of its four target genes, randomly selected among the genes supported by our transcriptome data. These novel qRT-PCR validated miR-128 target genes include EDARADD (EDAR-associated death domain), CAV1 (Caveolin 1) and DTYMK (deoxythymidylate kinase) (Figures 4C). One transcriptome-supported miR-128 target gene, RBM15B, could not be validated using qRT-PCR, as we did not observed significant changes in its expression level following miRNA precursor tranfection (data not shown).

Similarly, overexpressing miR-30 led to the downregulation of three of its five randomly selected transcriptome-supported targets. These qRT-PCR validated miR-30 targets include PLXNB (plexin B2), C11orf45 (chromosome 11 open reading frame 45) and ATP2B1 (ATPase, Ca++ transporting, plasma membrane 1) (Figures 4D). Two remaining miR-30 target genes, PPIA and CASD1 could not be validated via qRT-PCR, as we did not observed significant changes in its expression level following miRNA precursor transfection (data not shown).

Overall, at least 50% of the supported target genes tested in this study appeared to be regulated by their corresponding microRNA either directly or indirectly which is significantly better than the result of ordinary bioinformatic prediction with only 10% of correctly predicted target genes.

Functional analysis of MR-MIRNA target genes supported by transcrptome data

After functional classification of genes differentially expressed during myogenic differentiation in vitro, we thought to determine the proportion of microRNA-controlled genes within each functional class. Almost 30% of all protein coding genes contain miRNA seeds in their 3^′UTRs and are, therefore, potential miRNA targets [48]. Indeed, we have found that a significant proportion of genes differentially expressed during myogenic differentiation appeared to be controlled by microRNAs. However, the proportion of genes targeted by microRNAs was not the same within different functional classes. For example, 50% to 55% of genes related to angiogenesis and NF-kB signaling were potentially targeted by miRNAs, while in the case of myogenesis, potential miRNA targets represented only 8.6% of genes (Table 2). Interestingly, more miRNA targets were found among downregulated than upregulated genes. For example, among apoptosis-related genes, 17.6% of upregulated and 38.4% of downregulated genes were found to be targeted by miRNAs (Table 3). This observation prompted us to search for known functions of each supported target gene of every MR-miRNA described in this study using DAVID functional annotation tool [37, 38].

miR-1 and miR-206

It has been previously demonstrated that miR-1 and -206 target genes linked to the regulation of chromatin modifications [8], transcription [20, 49] and cell cycle [50]. Functional analysis of predicted miR-1 target genes supported by transcriptome data indicated that these microRNAs might also target genes involved in the regulation of apoptosis, DNA damage response, cell motility and protein modification, cell signaling and kinase activity [33]. Our own functional analysis of predicted target genes of miR-1 and -206 supported by transcriptome analysis confirmed that miR-1 and miR-206 might target genes involved in in these biological processes. In addition, miRNA-1 and -206 might target genes implicated in cytoskeleton organization, ubiquitination and nucleotide biosynthesis (Figure 5, Additional file 9).

Cell culture conditions and transfection

Primary human myoblasts were isolated from skeletal muscles of healthy subjects as described in [35], for details see Additional file 4: Table S4), purified with an immuno-magnetic sorting system (Miltenyi Biotec, USA) using an anti-CD56/NCAM antibody according to the manufacturer’s instructions. CD56-positive myoblasts were seeded in collagen-coated Petri dishes (P1) and cultured in DMEM, 10% FCS, 1% Ultroser G, at 37°C with 5% CO₂.All experiments were carried out between P1 and P5 to avoid cell senescence. Myoblast purity was determined by staining for Desmin. Purified myoblasts were plated in collagen-coated Petri dishes and cultured in growth medium containing DMEM supplemented with 20% fetal bovine serum at 37°C in humidified atmosphere with 5% CO₂. Myogenic differentiation of confluent cells was induced after 5 days by changing to DMEM containing 2% FBS (differentiation medium). Cells were kept in differentiation medium for 3 days. Myogenic fusion index at day three was 37-70%. Human immortalized myoblasts (iMyo) (kind gift of Dr. V. Mouly) were grown and differentiated as described in [113]. RD and TE671 (a kind gift of Dr. S. Leibowitz) were grown as described [114]. phrGFP-1 vector was from Stratagen, pcDNA3.2/V5 hsa-mir-128 (#26308) [115] and pCMV-miR30 (#20875) [116] vectors were provided by Addgene. Transient transfection of RD cells was performed in 6-well plate format using Lipofectamine 2000 according to manufacturer’s instructions with minor modification: 600.000 cells were added to the transfection mixture prepared directly in the cell culture plate. To transfect human immortallized myoblasts, anti-sense LNA inhibitors of miR-1, -206, 133a and b, and the irrelevant control were purchased from Exiqon (Denmark) and were mixed with LipoRNAiMAX (Invitrogen) in 6-well plates at 100nM final concentration, using 500 000 cells per well.

Reverse transcription and real-time PCR

For microRNA expression analysis 400 ng of total RNA purified via Trizol (Invitrogen) was converted into cDNA using 8 independent pools of primers (#4384791, Applied Biosystems, AB) and TaqMan microRNA Reverse transcription kit (#4366596, AB). cDNA was quantified using via qPCR using TaqMan 2x Universal PCR Master Mix, No AmpErase UNG (#4324018, AB) and human microRNA panel version 1.0 TLDA (TaqMan Low Density Array, AB), data were acquired on AB7900HT Real-Time PCR machine. Sequences of microRNAs and corresponding codes of TaqMan assays can be found in Additional file 2. For mRNA expression analysis total RNA was isolated from 2x10⁶ proliferating myoblasts, or differentiated myotubes using Trizol (Invitrogen) and reverse transcribed using the High Capacity cDNA Archive kit (Applied Biosystems, AB) according to the manufacturer protocol. cDNA was mixed with 2x Taqman PCR mix (AB) and amplified customized TLDA (TaqMan Low Density Array) using Abiprism 7900HT apparatus (AB). In case of miRNA and mRNA the expression level was calculated from Ct values using ΔΔCt method using GAPDH as a control [117]. The calculation of SEM in case of average of multiple samples has been done as described in [118].

Transcriptome profiling and data processing

Human primary myoblasts were sacrificed directly on plates at 30% confluency using Trizol, RNA was prepared using organic extraction and ethanol precipitation as described [119] followed by silica column cleanup on silica columns (Nucleospin RNA Extraction kit, Macherey Nagel). RNA extracted from individual myoblast lines was Cy5-labeled, mixed with with a pool of RNA samples labeled with Cy3 and hybridized to Gene Expression microarrays (4x44k #G4112F Agilent) and scanned as instructed by the manufacturer. Scanned images were then analysed using the Feature Extraction Software (Agilent), the subsequent gene expression data treatment was performed using R and Bioconductor. Spots with intensity lower than 50 or lower than background in more than 50% of biological replicates have been removed from further analysis. The background correction and intensity normalization procedures were applied for the remaining ~30 000 probes using the Bioconductor package vsn [120]. A background offset and a scaling factor for each array and dye channel were calculated using the least squares regression procedure, then the generalized log-transformation was applied. The ordinary least squares regression approach is based on the assumption that “most genes are not differentially expressed”. However, in the case of myogenic differentiation, where many genes are differentially expressed, this assumption does not hold. Therefore, to apply the above mentioned approach to myogenic differentiation, the vsn least squares model has first been applied to a subset of features, and then extended for the whole set of features. To select the subset of features, a pool of samples prepared from five different proliferating myoblast lines before and after myogenic differentiation were hybridized to two additional microarrays. Then 14.358 features that did not exceed the cutoff value of 1.23 fold change between the pools have been selected. To determine the differentially expressed genes, a t-test analysis was conducted using the limma package from Bioconductor [120]. Using this package, a linear model is fitted to the expression data for each gene. An empirical Bayes moderation of the standard errors were performed. This method borrows information across genes in order to arrive at more stable estimates of each individual gene’s variance, even for experiments with small number of arrays. The microarray data related to this paper have been submitted to the Array Express data repository at the European Bioinformatics Institute (http://www.ebi.ac.uk/arrayexpress) under the accession number E-MTAB-1577 (release date 2014-3-22).

miRNA targets predictions

RNA22 [40] and TargetScan algorithms [121] have been used to predict miRNA target genes in this study. miRNA targets identified by RNA22 have been downloaded from the RNA22 website http://cbcsrv.watson.ibm.com/rna22.html. Predictions were processed in the form of 3^′UTR binding sites occurrence matrices, containing the number of 3^′UTR binding sites for each combination of miRNA (columns) and mRNA (rows).

Pearson correlation between mRNA and miRNA expression profiles and permutation tests for significance

All the following has been performed using the R software. Pearson correlation values were computed using expression levels of each single transcript available and for each miRNA, found to be differentially expressed during myogenic differentiation in this study. Pearson correlation values were retained for further analysis, if they were corresponding to miRNA and mRNA pares having at least one occurrence in at least one of the two 3^′UTR binding sites occurrence matrices. To distinguish between biologically significant and random correlations, we performed a permutation test. We performed nine different permutations by calculating correlation values using a miRNA table with permuted sample designations. Then, the distribution of the real and permuted data were graphically plotted. For each miRNA, the 5% percentile of the Pearson correlation obtained with the permuted data was used as a threshold (th) for selecting biologically significant target genes. For each differentially expressed miRNA, the Pearson correlation value is compared to th and if smaller, the gene is retained as a biologically significant target gene regulated by the miRNA, with a 5% risk of error.