Identification of the miRNAome of early mesoderm progenitor cells and cardiomyocytes derived from human pluripotent stem cells

MicroRNAs are small non-coding RNAs involved in post-transcriptional regulation of gene expression related to many cellular functions. We performed a small-RNAseq analysis of cardiac differentiation from pluripotent stem cells. Our analyses identified some new aspects about microRNA expression in this differentiation process. First, we described a dynamic expression profile of microRNAs where some of them are clustered according to their expression level. Second, we described the extensive network of isomiRs and ADAR modifications. Third, we identified the microRNAs families and clusters involved in the establishment of cardiac lineage and define the mirRNAome based on these groups. Finally, we were able to determine a more accurate miRNAome associated with cardiomyocytes by comparing the expressed microRNAs with other mature cells. MicroRNAs exert their effect in a complex and interconnected way, making necessary a global analysis to better understand their role. Our data expands the knowledge of microRNAs and their implications in cardiomyogenesis.

cleavage sites of Drosha or Dicer in the microRNA maturation process, up to RNA editing by nucleotidyl transferases and exoribonucleases 12 . Variations could be present in both the microRNA 3′ or 5′-ends, although the latter are less frequent due to the possible changes in the seed sequence, and hence, in their complementary targets. However, 3′-end alterations have also been implicated in microRNA functionality, as they are thought to alter their biostability. The relevance of isomiRs is still debated but there is increasing evidence that they might have an important role in cell identity and behavior as they have specific expression patterns which may lead to a bigger impact in the mRNA regulatory network 12,13 . It has been proposed that the variations in sequence composition and length not only affect the mRNA target repertoire of the microRNAs, but also their stability and lifespan. Moreover, the isomiRs could interact with their homologue microRNAs, modulating their function and leading to a more efficient regulatory pattern 14 .
The potent ability of microRNAs in cell function regulation is given by a complex network formed by each individual microRNA. Potency is given by a double mechanism: on the one hand, a microRNA can target many mRNAs; whereas on the other hand, the 3′UTR of a mRNA can be targeted by many microRNAs. The whole mass of microRNAs then form a complex network that has been called the miRNAome of a cell, where there is a current set of microRNAs articulated in a given time to establish a specific cell phenotype. Small RNA-seq allows a detailed analysis of the microRNAs present at a given time. However, proper interpretation of these results require a deep analysis taking into account all the microRNAs and their variants, such as the isomiRs. With this perspective we analyzed the miRNAome of cardiac differentiation from pluripotent stem cells. We aimed to identify microRNAs networks involved in the regulation of this differentiation process. Instead of stressing the results of unique individual microRNAs, we proposed an analysis based on microRNA genome clusters and families.

RNA-Seq results: A general miRNAome at each stage of differentiation.
MicroRNA expression was analyzed by small RNA-sequencing in three biological replicates of three cell populations: pluripotent stem cells (PSC), primitive differentiated mesoendoderm progenitor cells (MPC 15 ), and cardiomyocytes (CM) (Supplemental Fig. 1). On average, approximately 42% of the total sequences were mapped to microRNAs, while the remaining reads were related to other non-coding RNAs, including rRNA, snRNA, snoRNA, and lincRNA (Supplemental Fig. 2). We found 694 microRNAs expressed in at least one of these cell populations (arbitrarily based on the expression of at least 5 reads in each replicate, and an average of 5 or more in at least one cell group). This is a significant part of the approximately 2800 human mature microRNAs described in mirBase database ( Fig. 1A; Supplemental File 1 (all microRNAs) and 2 (expressed microRNAs)). The most highly expressed microRNAs in PSC and CM are shown in Tables 1 and 2 respectively. The microRNA expression pattern significantly varied across cell populations. Although intra-group expression levels certainly showed some variations from sample to sample, a consistent expression pattern can be observed in the heatmap where some microRNAs increased or decreased throughout differentiation. We sought to further characterize this dynamic behavior during the three stages of differentiation by performing a fuzzy hierarchical analysis 16,17 . Fuzzy plots obtained by this method can be seen in Fig. 1B. After optimizing cluster adjusting (Supplemental Fig. 3) we found five different expression patterns of microRNAs. First, some microRNAs increased early once the differentiation had initiated and continued increasing up to the later cardiac stage. Second, some microRNAs increased at the later stage. Third, some other microRNAs decreased at the end of differentiation. Fourth, some rapidly decreased once differentiation was initiated. Finally, a group of microRNAs increased during the MPC stage, but decreased at the cardiac stage. MicroRNAs represented in each of these clusters can be found in the Supplemental File 3.
We performed a differential expression analysis between each cell population. As expected, many microRNAs were differentially expressed in each cell group. Volcano plots were generated by comparing the three cell populations in pairwise ( Fig. 1C; Supplemental Files 4-6). CM population showed the largest microRNA expression difference when comparing both PSC and MPC. There was a greater overlapping in the microRNA repertoire between PSC and MPC as they are closer in the developmental time course. For example, mir-302 family was highly expressed in PSC but maintained an upregulated expression at the mesoderm stage, with a significant decrease thereafter. These differences among groups were also seen in a distance heatmap (Fig. 1D). In agreement with the volcano plots, differences were more pronounced at the later stage of differentiation (cardiac stage) than between PSC and MPC. Finally, a principal component analysis (PCA) showed a marked independence of each cell group (Fig. 1E), stressing the significant differences among cell populations and validating the replicates. isomiR analysis. miRDeep2 software does not report variants from the published canonical microRNA sequences of mirBase, although it provides a pdf file for each microRNA detailing the aligned sequences (see Supplemental File 7 as an example). Clearly, many reads present a sequence that is slightly different from the exact microRNA. However, any sequence with more than one mismatch or a sequence with a larger extension than two bases in the 5′-end and more than five bases in the 3′-end are discharged from the analysis, and all the sequences that fulfill with the parameters are recorded in the alignment phase as part of the same microRNA molecule. There is accumulating evidence that suggest that these minor modifications in the 5′-and 3′-arms may be functional, and more importantly, their gene targets may or may not overlap with the canonical microRNA targets, hence it has been suggested that a proper analysis of the miRNAome should include the study of isomiRs 18 . In the last years several packages have been developed for the analysis of isomiRs. We used the miraligner/isomiR package to detect these isoforms in our samples as it has been shown to be highly accurate 19,20 , although we also ran the analysis using Chimira and Isomir-SEA 21,22 . miraligner/isomiRs packages searches for matchings of the exact microRNA sequences and distinguished those with modifications in the 5′-end and in the 3′-end (up to 3 extensions or deletions in the extremes). Figure 2A shows the non-filtered count of exact sequences and those with modifications in either one or both ends. Overall, 6,359 sequences were detected (Supplemental File 8). After filtering by a mean of 5 or more in at least one cell group (as in the previous analysis), there were 2809 sequences. Modifications were frequent, and even the exact sequence was in fact less numerous than the modified one. This finding was consistent when the analysis was run with other softwares, with even more modifications beyond position 3 (Supplemental Fig. 4).
We then examined the expression levels of the isomiRs in comparison to exact microRNAs. Figure 2B shows that many isomiRs were highly expressed, even to the level of the exact microRNAs. There were no differences between the expression levels among cell groups (p = 0.18 by ANOVA). Figure 2C shows the position of the isomiRs modification along the microRNA sequences. As it can be seen, there were no differences among cell groups. When considering modifications at the 5p-end, most changes occurred at position -3 (that is, deletion of 3 bases). In the case of the 3p-end, most modifications occurred at position -1. In any case, modifications were frequent at all positions. Nucleotide substitutions ("SNPs") occurred at a low percentage. However, it is interesting to note that substitutions were more frequent at the 5p-end extreme (precisely at the seed) (Fig. 2D). Moreover, these nucleotide substitutions involved a particular subset of them, including C > G, G > A and T > C (Fig. 2E). Both findings did not occur at random, and hence, we believe that they correspond to real base substitution and not sequencing errors.
The highly expressed isomiRs could have an impact on the regulation of gene expression producing an important transcriptomic effect, particularly those with variations in the 5′-end. As an example, we described an interesting (and extreme) case of a seed-modificated isomiR. mir-302 family members are highly expressed microRNAs in PSC. We noticed a higher expression of several isomiRs related to hsa-mir-302a, particularly at the 5′-end. Figure 2F shows the alignment of isomiRs and exact sequences, as well as the expression level. As can be noted, the seed is significantly changed in these isomiRs. The potential targets of these sequences were analyzed (Fig. 2G) by prediction using the DIANA tools (MR-microT software 23 ,). Interestingly, for the miR-302a-3p there were several common targets between one of the isomiRs and the microRNA exact sequence, but many were exclusive for each isoform (94 for 302a3p-iso1, 81 for 302a3p-iso2 and 181 for 302a3p-iso3). Similarly, isoforms of the miR-302a-5p had several unique targets (242) that were not present in the repertoire of targets regulated by the exact microRNA. To explore whether the isomiRs regulate genes from different biological processes, we analyzed the gene ontology terms (GO) associated with targets specific to the exact microRNA and targets common to both microRNA and isomiRs. Surprisingly, the analysis of the GO terms including all predicted targets (from both microRNA and isomiRs), revealed an enrichment in new biological processes not found by the analysis of the exact microRNA targets alone (p < 0.01), such as the Wnt signaling pathway, calcium modulating pathway, or regulation of cellular component size (Supplemental Fig. 5). Hence, isomiRs potentially expand and strengthen the regulatory network of the microRNAs. Finally, we evaluated the correlation between the target affinity scores provided by MR-microT of the exact microRNA and isomiRs. A linearity in this regression means that there is a correlation between the probability of a gene to be regulated by both microRNA and respective isomiR. Interestingly, we observed that none of the analyzed pairs showed a high correlation, implying that even though the microRNA and isomiRs might have some overlap in the regulated genes they present different target affinities,  hence denoting that they potentially exert their inhibitory effects in an additive and complementary manner more than in a redundant way (Supplemental Fig. 6). MicroRNA modifications also occur by a different mechanism. A-to-I edition or ADAR (adenosine deaminase acting on RNA) are identified in short RNAs. These modifications are believed to alter microRNA degradation 24,25 . We used the Chimira software from the EBI institute 21 to identify ADAR modifications and found consistent alterations in the microRNA sequence. ADAR modifications were globally found across all microRNA hairpins, in a similar fashion as isomiRs, on each extreme as well as in the internal structure of the mature microRNA (Fig. 3A). Proportionally, we found a significantly greater number of ADAR modifications in the PSC population compared to CM (p = 0.014, Fig. 3B). When we analyzed the expression of the ADAR-modified microRNAs, they were similarly expressed across differentiation (Fig. 3C). However, the expression level was low, as it has been described previously for ADARs (Fig. 3D) 25 . Altogether, these results suggest that isomiRs may have a functional significance in pluripotent stem cells and during the cardiac differentiation. Their equivalent number of reads to the exact sequence microRNAs and their different targeting specificities support this concept. With the study and incorporation of the isomiRs in the analysis, we determined a miRNAome during cardiac differentiation far richer and more complex than previously thought. microRNA genomic cluster and family analyses. Even though unique microRNAs are associated with specific cell functions, it is increasingly clear that miRNAome acts in concert. Non-physiological experiments may mask the real role of each microRNA in the differentiation process, such as a forced increased expression of a single microRNA to supra-physiological levels. We think a better understanding of the function of the miR-NAome is then given by analyzing groups of microRNAs as a whole. microRNAs are usually grouped based on two features: genomic clustering, based on the genomic proximity (arbitrarily defined in mirBase as those microRNAs located within 10,000 bases), and seed sequence, based on sequence similarity. In many cases the genomic distribution of the microRNAs is not random and it has been found that some of them group in clusters in specific genomic regions and aggregate in islands of microRNAs. It is likely that these islands share regulatory machineries, although little is known about them. However, all the microRNAs that share the same seed sequence constitute a microRNA family, and are expected to regulate similar gene targets. We therefore analyzed the expression profiles of the microRNAs during cardiac differentiation considering clusters and families.
A first insight of these analysis is shown in Fig. 4A. The Circos plot shows the genomic distribution of the microRNA clusters. Each ring contains bar graphs corresponding to the average expression of the microRNAs clusters of PSC (inner circle), MPC (middle circle) and CM (outer circle). In the center region, the linking lines represent the association between clusters and families. As it can be seen, many clusters are composed of microR-NAs from different families (blue crossing-lines in the center). For example, four families (mir-17, mir-19, mir-25 and mir-363) are located at three different clusters (at cluster 7.4, cluster 13.2, and cluster X.8). Many other clusters are composed by a single family (red lines), as for example the mir-302 family/cluster (cluster 4.1).
By comparing heatmaps of both cluster and family groups, it can be seen that they were quite similar (Fig. 4B). This finding highlights the fact that clusters and families share a similar dynamic behavior during differentiation, even though they may be spread along the complete genome. Expression level for clusters and families are shown in the Supplemental Files 9 and 10. The analysis of similarities between cell populations by a sample distance heatmap (Supplemental Fig. 7A), showed a high similarity between the cluster and family. We also performed principal component analysis on the microRNA profiles to assess the differences and similarities between the three cell populations and these results support the concept that microRNA families and clusters are very similar, therefore being more adequate to discriminate between cell populations than individual microRNAs (Supplemental Fig. 7B).
We then studied the predicted targets by gene ontology of three cluster/families of microRNAs specific to each stage of the differentiation process to evaluate the pathways they might be regulating during cardiac differentiation. We chose the miR-302 family, miR-17/92 family and a group of microRNAs highly expressed in CM population. Figure 5A shows the number of predicted targets that are regulated by one or more microRNAs from each group. Many genes are under the control of more than one microRNA, which illustrates the cooperative role that microRNAs have in gene silencing. We next analyzed the Gene Ontology (GO) terms associated to the predicted target genes to study the putative molecular processes regulated by these groups of microRNAs (Fig. 5B). For both miR-302 and mir-17-92a families many GO terms were related to developmental processes (see complete lists in Supplemental Files 11 and 12). It is worth to note that heart development ranked at the top in both families. As both microRNAs families were highly expressed in PSC, it was expected that these regulated GO pathways were expected. When microRNAs where individually analyzed, however, GO terms did not clearly relate to embryonic development (see an example in Supplemental Fig. 8), emphasizing that the addition of all family members predicts a highly significant effect over embryonic development, particularly cardiac differentiation. We then analyzed the GO terms associated to targets of the microRNA families highly expressed in CM, including miR-1 (mir-1 and mir-206), miR-133 (miR-133a and b), miR-208 (miR-208a and b), miR-490, miR-499 and miR-143. The analysis revealed that these microRNAs might be downregulating genes related to different cell lineages derived from the mesoderm germ layer, such as striated muscle development, skeletal system development and urogenital system development (see complete list in Supplemental File 13). This is in line with the concept that microRNAs act as key regulators during early development and orchestrate the gene expression that leads to the establishment of cell identity, and in the specific case of CM, inhibiting the development of other mesoderm cells.
Additionally, we studied the correlation between the target affinity scores of the miR-302 cluster/family for their common targets (see above). In Fig. 5C the upper part of the matrix shows the correlation scores between pairwise comparisons of miR-302 members. As it can be observed, -3p microRNAs had the highest correlation between them than -5p microRNAs, except for miR-302b-5p and miR-302d-5p. This observation implies that Scientific RepoRts | (2018) 8:8072 | DOI:10.1038/s41598-018-26156-3 members of this cluster act cooperatively in the control of gene expression, but redundancy between them is not equal for all the members of the cluster. This latter observation suggests that despite the overlapping of gene targets, each microRNA from miR-302 cluster interacts with different affinities to their target RNAs. Similar findings for the mir-17-92a cluster/family and the cardiac microRNAs can be seen in Supplemental Fig. 9. Hence, these findings support the previous concept that embryonic microRNA cluster/families team up to inhibit different developmental pathways.
Cardiac associated microRNAs. Comparing microRNAs expressed in cardiac cells to pluripotent stem cells does not confirm that those differentially expressed microRNAs are exclusively of cardiac cells. Hence, we compared the miRNAome of cardiac cells against other differentiated cells. An external small-RNA-Seq database was analyzed with the same methods (GEO accession number GSE68189) 26 . These samples were mostly from the same embryonic stem cell background (H9), and includes mesenchymal stem cells derived from PSC and  three neural differentiated cells (fore-, mid-, and hind-brain). This analysis also included a cell line derived from human embryonic lung fibroblasts. There was a correlation of the miRNAome between PSC samples, hindered by some microRNAs with low expression in either group (Fig. 6A). A heatmap of the microRNA expression profiles of all cell populations showed a marked difference between the differentiated cells (Fig. 6B). As it can be seen, a microRNA expression pattern was distinctive. Pairwise comparisons could be observed in the volcano plots (Fig. 6C). Supplemental Files 14-18 list the microRNAs down-regulated and up-regulated in each comparison between differentiated cells. Finally, we analyzed the expression of microRNAs comparing all cell lines simultaneously in order to build a cardiac miRNAome. In Table 3 we show three different microRNA groups. The first two columns show the top ranking microRNAs with highest absolute expression in CM. The third and fourth columns show that those microRNAs are significantly up regulated in CM compared to all other cell lines analyzed, including PSC and differentiated cells. Finally, the last two columns show those microRNAs uniquely expressed in CM (that is, zero expression was found in any other cell population).

Discussion
Cardiac differentiation from human pluripotent stem cells presents dynamic changes in the landscape of microRNA expression. We analyzed the microRNA signature in order to establish the identity of pluripotent stem cells (PSC), early mesoderm progenitor cells (MPC) and cardiomyocytes by small RNA-sequencing. Our analysis revealed that almost 700 microRNAs are differentially expressed in the three cell populations and they have specific expression patterns according to the differentiation stage. Therefore, each cell population has a characteristic miRNAome that has an important regulatory role in the establishment of cardiac lineage. Two additional concepts were explored in this paper. First, families and clusters of microRNAs work in harmony through the dynamic stages of differentiation. Second, minor but numerous isoforms of the microRNAs may take a significant role in the regulation of this differentiation.
MicroRNAs form a complex regulatory network important in the control of gene expression during human embryonic development. There is increasing evidence of the role they have in the post-transcriptional regulation of key developmental genes 27,28 . The complexity that characterizes the way microRNAs work is partially explained by the fact that one mRNA is simultaneously targeted by numerous microRNAs and each microRNA regulates several mRNA at the same time. Moreover, the high sensitivity in new sequencing technologies allowed the detection of microRNA isoforms named isomiRs. These variations imply an extra level of regulatory complexity as each microRNA can have a myriad of isomiRs that might strengthen and expand the regulation of the target RNA messengers. The isomiRs were initially considered errors from either the sequencers or softwares used in the analysis of the results 12 , whereas today there is evidence that supports their functional relevance 29,30 . We described an extensive network of isomiRs during cardiac differentiation. Most abundant isomiRs were those with variations in the 3′-end. Nonetheless, we also identified several isomiRs with changes in the sequence at the 5′-end. These modifications in the 5′-end presumably have an impact in mRNA recognition. Despite this, it has been proposed that heterogeneity in the 3′-end might have important implications in the microRNA stability and AGO loading affinity 12 , hence having a potential effect in the microRNA regulatory activity. Importantly, we found no evidence that differentiation affects isomiRs formation from the exact sequences. We then analyzed in more detail the microRNA miR-302a and their isomiRs, revealing the case of a particular microRNA being less abundant than its isomiRs and with a repertoire of predicted targets considerably different for most of the discovered 5′-isomiRs. We identified 415 novel predicted targets for the 5′-isomiRs of the miR-302a-3p and 242 for the 5′-isomiRs of the miR-302a-5p, both highly expressed in the PSC. Furthermore, the Gene Ontology analysis of common and different transcripts targeted uncovered new cellular processes that are implicated principally in cell differentiation, which would not have been identified if the 5′-isomiRs were not considered in the analysis. Our data is in line with the work 13 that demonstrates that the incorporation of the isomiRs increases the accuracy of target prediction.
MiRNAome analysis has generally overemphasized individual microRNAs expression profiles. However, the sequence of many microRNAs are very similar, particularly in the seed region, which has been used to group them into families. Many microRNAs are also typically found in genomic clusters, arbitrarily defined as microR-NAs separated by no more of 10,000 bases. We then believe that these features should be considered when analyzing a miRNAome, as families and clusters may exert their function over common genes or common molecular pathways. Therefore, we determined the dynamic expression profiles of 63 microRNA clusters and 141 microR-NAs families that contribute to cell fate decision during cardiac differentiation. The most expressed microRNA family in PSC and MPC was miR-302, which also corresponds to a cluster we named Chr4.1 (first one found in Chormosome 4). This family/cluster has been extensively studied in previous works 31,32 and it is of great importance in the maintenance and induction of pluripotency. In addition, there is another highly regulated cluster known as C19MC located at chromosome 19 (Chr19.10). This cluster is one of the largest found cluster in primates and has been suggested in a recent work to regulate apoptosis in pluripotent stem cells 33 . This cluster includes the family miR-515 but also miR-373 and miR-290. Both families are also up regulated in PSC. Another  interesting upregulated cluster in both PSC and MPC is located on chromosome 13 (Chr13.2), which includes members of the families miR-17, miR-19 and miR-25. These microRNAs have an important role in normal tissue development but also an oncogene role in different cell types 34,35 . Another family expressed in MPC that stands out from the rest is the miR-26 family (miR-26a and b), which is also related to tumorigenesis 36 . These groups of microRNAs could be associated with the epithelial to mesenchymal transition occurring in the development of the mesoderm progenitor cell, which is found in developing tumors. Lastly, we looked at the up-regulated microRNA clusters and families in CM. The five clusters with higher expression level in CM are Chr5.2, Chr11.7, Chr19.8, Chr8.7 and Chr1.2. By large, these clusters are formed by members of the miR-143 family, let-7 family and miR-30 family. Consistently with previous works, miR-143 regulates smooth muscle differentiation and let-7 family decreases with early differentiation and is up regulated in adult cardiac tissue 37,38 . Gene ontology analysis of the clusters/families in the cell groups showed that there are shared processes regulated such as heart development or muscle development, and most GO terms are related to embryonic development. Notably, there is a progressive narrowing (in terms of embryo development) on the range of cellular processes regulated once the cells have reached a more mature stage. In PSC there are more GO terms related to a wider variety of cell lineages like heart, skeletal system, muscle tissue, forebrain, striated muscle tissue, in utero embryonic development, mesenchymal, telencephalon and pallium development. On the contrary, GO terms related to genes silenced in MPC and CM are mostly general cellular processes, neural and muscle development. This observation suggests that once the cells have reached a mature phenotype,as in the case of a cardiac phenotype, the up-regulated microRNAs control gene expression of very specific pathways. Moreover, the score correlation of predicted targets of the family miR-302 illustrates the synergistic effect that these microRNAs have. Finally, it is remarkable that the individual GO analysis of a microRNA does not usually show enough gene targets of embryo development. It is necessary to include all family members (as in the case of mir-302) in order to reach significance in the developmental pathways. Altogether, these results suggest that microRNAs act cooperatively in silencing key genes during cardiac differentiation, emphasizing the importance of a genome-wide analysis, in addition to considering the expression profile of families and clusters of microRNAs.
Finally, we compared our results with published data of microRNAs in different cell types to narrow the list of microRNAs that specifically contribute to the miRNAome of the CM (Box 1). Our analysis showed that there are 848 microRNAs differentially expressed between the six cell types. In line with previous observations, each cell type has a specific signature of microRNAs implying the importance of microRNAs in cell fate decision. Lastly, in our analysis we were able to identify a group of microRNAs that are significantly up regulated in cardiac cells and a short series of microRNAs that are unique to CM. This finding is relevant, as they confer CM with a specific signature. Some of these microRNAS have been previously reported and are well-known in the context of cardiac differentiation, such as hsa-mir-1-3p, hsa-mir-499-5p, and the mir-208 and mir-133 families 39,40 . However, other microRNAs have not been described in cardiac differentiation, including hsa-mir-675-5p, hsa-mir-585-3p, hsa-mir-1248 and hsa-mir-509-5p. To what extent these microRNAs are relevant in cardiac differentiation still remains to be determined.
In conclusion, we determined the microRNAs that are involved in cell lineage specification during human pluripotent stem cells differentiation towards mesoderm and cardiac cells. These results expand the knowledge of microRNAs and their isomiRs in the regulation of genes and molecular mechanisms important during cell fate decision. Due to the high complexity of the miRNAome, a genome-wide approach is mandatory for a better comprehension of the role these small non-coding RNAs have and functional studies that include all microRNAs members of families/clusters will be necessary to further understand their proper relevance in the differentiation process.

Methods
Cell Culture and cell differentiation. Human embryonic stem cells (H9-hTnnTZ-pGZ-D2) were obtained from WiCell and maintained in co-culture with irradiated primary mouse embryonic fibroblasts (MEFs) in the presence of FGF (8 ng/ml; Thermo Fisher Scientific) under conditions described by the supplier, as previously reported 41 . This cell line has been genetically modified in WiCell to express both GFP and Zeocin resistance under the promoter of cardiac troponin T. All experiments were carried out in accordance with relevant guidelines and regulations. All cultures were maintained at 37 °C in a saturated atmosphere of 95% air and 5% CO2. To induce mesoderm differentiation and obtain the early mesoderm progenitor (MPC) we followed a previously described protocol by Evseenko et al. (Supplementary Fig. S1A) 15 . Briefly, cells grown in iMEFs were detached using TrypLE TM Select (1X; Thermo Fisher Scientific) and cultured on 6-well plates pre coated with Geltrex Matrix (Diluted 1:1000 from 15 mg/ml; Thermo Fisher Scientific) with mTesR (StemCells Technologies). Cells were grown in these conditions until they were 60-80% confluent (typically 3 days). To induce differentiation medium was replaced with StemPro ® -34 SFM (Thermo Fisher Scientific) supplemented with the following proteins: human bone morphogenic protein 4 (BMP4), human vascular endothelial growth factor (VEGF) and fibroblast growth factor 2 (bFGF) (all at 10 ng/ml; Thermo Fisher Scientific) with the inclusion of Activin A only at the first day (10 ng/ml). At day 3.5 after mesoderm induction, the MPC were isolated by fluorescence-activated cell sorting (FACS) using PerCP/Cy5.5 anti-human CD326 and APC anti-human CD56 antibodies from Biolegend. For CM differentiation, embryoid bodies (EB) formation was forced by centrifugation of hESC (H9-hTnnTZ-pGZ-D2) in 96-well plates V-bottom with mTeseR supplemented with the growth factor BMP-4 (10 ng/ml; Thermo Fisher Scientific). After 24hs, EB were cultured in suspension with StemPro ® -34 SFM (Thermo Fisher Scientific) supplemented with BMP-4 (10 ng/ml), bFGF (5 ng/ml) and Activin A (3 ng/ml) during 4 days. At day 4, EBs were collected and passed to adherent plates precoated with gelatin (0,1%; Sigma Aldrich) and the inhibitor Wnt response factor (IWR-1) was added to the medium (5 microM; Sigma-Aldrich). Between day 5 and day 21 cells were cultured in StemPro ® -34 SFM supplemented with VEGF (5 ng/ml), bFGF (10 ng/ml) and IWR-1 (5 microM). At day 11 CMs were purified with two pulses of Zeocin for 72hs (75 microg/ml), allowing 2-3 days to recover between pulses. At day 21 the percentage of cells GFP+ was analyzed by Flow Cytometry. At this stage, cardiomyocytes are differentiated up to the point of presenting spontaneous beating activity.
Small RNA-Sequencing samples preparation. The small RNA molecules were extracted with mirVana TM miRNA Isolation Kit (Ambion, Thermo Fisher Scientific) following manufacturer's protocols. Small RNA libraries were prepared with 200 ng of RNA and using the NEBNext Small RNA Library Prep Set with modified adaptors and primers compatible for Illumina (New England Biolabs) and cleanup and size selection (50pb) was made by gel electrophoresis following manufacturer's protocols. Sequencing was carried out at the TCGB Resources -UCLA Path & Lab Med using an Illumina HiSeq. 2500, in one lane and single read sequencing. Sequencing datasets are deposited in GEO (GSE108021).
Small RNA-sequencing results analysis. Processing and annotation of sequences based on identity to known transcribed RNAs or as novel miRNAs was performed by Novogenix BSU using the CAP-miRSeq (v1.1) pipeline 42 . Briefly, pretrimmed and posttrimmed reads were screened for quality, specificity of mapping, and contaminant sequences using FASTQC. Then prior to alignment, low-quality bases, homopolymer sequences, sequences matching the first 3 bp, and the reverse complement of the adaptor sequences for Illumina were trimmed using Cutadapt (v1.13) with a maximum error rate of 0.1%. The trimmed reads were then aligned using Bowtie and HTSeq 43,44 for the reads counting and then the reads were mapped to the human genome reference (GRCh38) with miRDeep2 mapper and to miRbase (v21) with miRDeep2 module which allowed us to identify known and novel microRNAs 45 . Normalization of read number was carried out as Counts per Million within each sample. For differential expression analysis we used R package DESeq. 2, using raw data as inputs. IsomiR analysis was done mainly with miraligner/isomiR (v1.6) 20 , but also with the isomiR-SEA 22 and Chimira web interface 21 . Gene target prediction of microRNA was done with miRNAtap from Bioconductor package, which includes the prediction algorithms of DIANA, Miranda, PicTar, TargetScan and miRDB. Targets of isomiRs were predicted using MR-microT from DIANA Tools 46 . Gene Ontology analysis was done with enrichGO from ClusterProfiler package in R 47 .
Accession code. RNA-sequencing results can be found in GEO with the accession number: GSE108021.