{'fraction_non_alphanumeric': 0.033687710548190926, 'fraction_numerical': 0.022837642735267095, 'mean_word_length': 5.1166711265721165, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic regulation of small RNAome during the early stage of cardiac differentiation from pluripotent embryonic stem cells

Embryonic stem cells (mESCs), having potential to differentiate into three germ-layer cells including cardiomyocytes, shall be a perfect model to help understanding heart development. Here, using small RNA deep sequencing, we studied the small RNAome in the early stage of mouse cardiac differentiation. We found that the expression pattern of most microRNA (miRNA) were highly enriched at the beginning and declined thereafter, some were still insufficiently expressed on day 6, and most miRNAs recovered in the following days. When pluripotent embryonic stem cells are differentiating to cardiomyocytes, targeted genes are concentrated on TGF, WNT and cytoskeletal remodeling pathway. The pathway and network of dynamically changed target genes of the miRNAs at different time points were also investigated. Furthermore, we demonstrated that small rDNA-derived RNAs (srRNAs) were significantly up-regulated during differentiation, especially in stem cells. The pathways of srRNAs targeted genes were also presented. We described the existence and the differential expression of transfer RNA (tRNA), Piwi-interacting RNA (piRNA) and Endogenous siRNAs (endo-siRNAs) in this process. This study reports the genome-wide small RNAome profile, and provides a uniquely comprehensive view of the small RNA regulatory network that governs embryonic stem cell differentiation and cardiac development.


Introduction
Small RNAs are master regulators controlling proliferation and differentiation, particularly in regulating stem cell biology and cardiac development. The switch from pluripotent to lineage-specified cells is accompanied by the up-regulation of many small regulatory RNAs, such as microRNAs (miRNAs) [7] noncoding endogenous RNAs with approximately 18-25 nucleotides in length, which regulate gene expression mainly by inhibiting translation or promoting degradation of target mRNAs by base paring with specific mRNA targets [11].
Using combinations of next-generation sequencing (NGS) technologies, small RNA population was profiled in numerous areas of biological research, including health and diseases [10].
Lots of heart disease are associated with decreased functional cardiomyocytes. It is well established that embryonic stem cells (ESCs) are able to generate bona fide cardiomyocytes. ESCs derived from the inner cell mass (ICM) of preimplantation embryos can be propagated in vitro in an undifferentiated state and, when allowed to differentiate, can form endodermal, ectodermal, and mesodermal derivatives in vitro and in vivo. ESCs thus can be a potential source of donor cardiomyocytes (or alternatively, donor cardiomyogenic progenitors) for therapeutic interventions targeting heart diseases.
Cellular differentiation is a sequential commitment process achieved through multiple intermediate states, and is tightly regulated by different levels of gene expression program. Cardiac differentiation from ESCs represents a paradigm for studying cell fate determination in both cellular and molecular levels. Spontaneous cardiac differentiation from mESCs can mimic the in vivo cardiac development. Thus, understating how the cell fate of cardiac cells is determined not only holds unlimited potential for the cell replacement therapies, but also yields general principles of cell fate determination and organ development [22].
In vitro ESC differentiation typically requires an initial aggregation step to form structures, termed embryoid bodies (EBs), which further differentiate into cardiomyocytes [1]. Thus, in vitro cardiac differentiation involves several distinct processes, including stem cell aggregation to form embryoid bodies, stem cell commitment and final functional cardiogenesis, which entails maturation to form beating cardiomyocytes. It has been shown that miRNAs can promote cell cycle exit and differentiation, and are required for proper ESCs differentiation [7]. However, the mechanism through which small RNA contribute to stepwise differentiation is currently unknown, and a global view of small RNA repertoire is still lacking.
To comprehensively understand the dynamic changes of small RNA expression during cardiac differentiation, we applied next-generation sequencing to analyze the small RNA fraction on day 0 (D0), day 2 (D2), day 6 (D6) and day 9 (D9) of differentiation. These time points afforded a view of critical transcriptomic changes occurring both before and during cardiac differentiation from pluripotent stem cells.
Our results provide a global view of RNAome, including miRNAs, srRNAs, piRNAs, endo-siRNAs and tRNAs, as well as their expression pattern during differentiation. We also provide a potential resource for the regulatory role of miRNAs and srRNAs that target signaling pathway during the differentiation of ESCs to cardiomyocytes.
In vitro differentiation of mESCs to cardiomyocytes was performed using the hanging drop method [20]. Briefly, ESCs were trypsinized and resuspended in complete growth medium without LIF, and plated in hanging drops at around 400 cells. They were dropped per 20 μl on the lids of petri dishes filled with PBS and cultivated in hanging drops for two days. EBs then were collected and plated on low-attachment cell culture plates for two days [9]. They were further plated separately onto gelatin coated 24-microwell on D6. For the determination of cardiac differentiation, spontaneously beating cardiomyocytes were investigated from the outgrowth of EBs cultured in 24-microwell plates. The well harboring beating cells were marked and cells in it were picked up on D9. Cells were collected on D0, D2, D6 and D9.

Small RNA cloning and sequencing
To obtain high quality total RNA samples, embryoid bodies, spontaneously and rhythmically contracting cardiomyocytes were manually picked up under compound microscope. Total RNA from D0, D2, D6, and D9 were extracted using Trizol (Invitrogen). The quality and integrity of the total RNA was evaluated by electrophoresis on 1.2% agarose gel and Agilent 2100 BioAnalyzer (Agilent). Small RNAs ranging from 18 to 30 nt were gel-purified and ligated to the 3′ adaptor and 5′ adaptor oligonucleotides. Samples were allowed for deep sequencing on the Illumina GAII platform.

Bioinformatics analysis of small RNA data
In order to make statistics of miRNA reads, mirTools [23] was used. Reads were mapped to the University of California at Santa Cruz (UCSC) mm10 assembly of the murine genome by using megablast. miRNAs profiling was according to miRNA (miRBase version 18). Normalization formula is by using transcripts parts per million (TPM): normalized expression = actual miRNA count / total count of clean reads × 1,000,000. Here, unique reads refer to different types of reads, and redundant reads refer to total reads.
Mapped miRNA which has reads in at least three samples were further analyzed using Gene Cluster 3.0. Gene Tree view which can display hierarchical as well as k-means clustering results was used.
The miRNA targets prediction software was composed of validated and the computationally predicted targets. The validated target data sets were selected from TarBase, miRecords, whereas miRanda, TargetScan, PicTar, PITA, MirTarget2, RNAhybrid, RNA22 were used to predict the targets of miRNAs. The majority rule voting was selected to predict the targets. After that, the validated and the predicted targets are combined as miRNA targeted gene pool [11]. PACCMIT-CDS, RNAhybrid, miRanda and MicroTar were used in srRNAs targets prediction.

Functional annotation of target genes
To group miRNA targets that are functionally related into known functional categories and pathways, metaCore-GeneGo Pathway Maps were used [11]. For highly enriched miRNAs, their targets were predicted and uploaded to GeneGo. Targets were analyzed with the Pathway Maps [11]. Regulatory networks were drawn using Cytoscape (http://www.cytoscape.org).

Novel miRNA and endo-siRNA prediction
Small RNA sequences were mapped to the genome and fetch each exact sequence match along with 100 bases flanking either side. For predicting the secondary structure, mfold was utilized (http://mfold. rna.albany.edu/?q=mfold/RNA-Folding-Form).
Reads were aligned themselves by using megablast. We only fetched the plus/minus paired reads. Candidates were mapped to the transcriptome SRR1028908. By using mfold, the secondary structure of transcriptome was predicted. The paired reads so as to leave 3′, twonucleotide overhangs, and on the long double-stranded transcriptomes was identified as novel endo-siRNAs. To obtain a comprehensive overview of small RNA dynamics during cardiac differentiation, a feeder-free murine embryonic stem cells (mESCs)-derived cardiomyocyte differentiation program was used. Samples were collected from four time points: pluripotent embryonic stem cells (D0), undergoing EB formation (D2), early stage of differentiation (D6) and rhythmically contracting cardiomyocytes (D9). Cells were harvested individually followed by applying to the next generation sequencing platform for retrieving small RNA expression.
Summary of small RNA sequencing data is shown in Table S1A. Numbers of high quality reads at four time points are 11989936, 12984769, 13221096, and 13234931, which correspondingly represent 6545457, 5696149, 6775264, and 7227913 small RNA sequences. The percentage of small RNAs mapped to rfam (ribosomal fragments) varies significantly (7.76%, 74.93%, 46.55%, and 19.19%). Also, the percentages (44.98%, 7.8%, 16.37%, and 53.93%) of miRNAs changed dramatically during the differentiation. This data indicates that the miRNAs and rRNAs are two major small RNA populations in cardiac differentiation from embryonic stem cells.
Unique reads referring to different types of small RNA sequence were shown in Table S1B. Those corresponding to unannotated reads are relatively high (54.05%, 42.38%, 48.38%, and 50.84%), suggesting that there are still lots of unannotated small RNAs species remain to be determined.

Length distribution of small RNAs and miRNAs
Length distribution of total redundant small RNAs is shown in Fig.  S1A. Majority of small RNAs varies from 18 to 22 nucleotides in length, corresponding to the size of miRNAs. Length distribution of unique small RNAs is shown in Fig. S1B.
Given that miRNAs play a critical role in the direct stem cell differentiation [7], firstly we analyzed the character of known miRNAs. The length distribution of total redundant miRNAs is shown in Fig. S1C. It has a single peak at nucleotide 22, typical for miRNAs. However, these miRNAs do not account for the entire 22-nt peak in Fig. S1A, suggesting that a new population of small RNAs remains to be discovered. Length distribution of unique miRNAs is shown in Fig.  S1D. Interestingly, small RNA reads at D2 showed the lowest counts (Fig. S1), suggesting that early cardiac differentiation displays drastic changes in small RNA expression profile.

Expression profile of miRNAs during mESCs differentiation into cardiomyocytes
Next, we compared the changes of miRNA expression. All known miRNAs with normalized expression levels are shown in Table S2, while as 30 of the most abundant miRNAs in the four samples used in this work are presented in Table S3.
Small RNA-Seq generated a comprehensive and quantitative atlas of miRNA expression during ESCs differentiation into cardiac cells (Fig.  S2). Heatmap clustering indicates that miRNAs on D0 and D2 have a closer expression pattern, while the pattern of D6 and D9 are similar.
Subsets of miRNAs were mainly organized into four categories. Firstly, the miRNAs that were highly enriched on D0 and immediately declined during the following days (Fig. 1A). These miRNAs may have specific function in maintaining stem cell properties. Indeed, miR-290-295 family has been shown to have a critical role in maintaining pluripotent state by regulating self-renewal [8].
Secondly, 32 miRNAs were down-regulated on D2, but gradually recovered on D6 and D9 (Fig. 1B). The down-regulation of these miRNAs might be involved in de-repression of genes that initiate cell differentiation. In accordance with this speculation, it has previously been reported that miR-185 and miR-191 are involved in the differentiation of primary hippocampal neuron, and these miRNAs potentially play a role in the initial stages of the differentiation [18].
Thirdly, 14 miRNAs were down-regulated on both D2 and D6, and recovered on D9 (Fig. 1C). These miRNAs may have a function in promoting the appearance of the cardiac progenitor cells. Finally, 12 miRNAs up-regulated on D9 (Fig. 1D) and can thus potentially be biomarkers of mature cardiac cells. Besides, we also noted several other patterns, shown in Fig. S3.
Based on these analysis, the change pattern of the most miRNAs can be summarized: some highly enriched at the beginning and declined thereafter, some were still insufficiently expressed on D6, and most miRNAs recovered in the following days. The expression patterns of 9 miRNAs were verified with real-time quantitative PCR (Fig. S4) in mESCs (R1) and displayed similar expression patterns with the deep sequencing data.
Deep sequencing is especially effective in the discovery of novel transcripts. To discover novel miRNAs, sequences with at least 100 reads counts were selected. Through analyzing unannotated small RNAs with mfold, 25 novel miRNAs were identified. They have a typical hairpin loop secondary structure as shown in Fig. S5, which means there is still some low abundance of miRNAs to be discovered. Sequences of these novel miRNAs were shown in Table S4.

Pathway enrichment analysis of special miRNAs targets
To identify the possible biological functions of differentially expressed miRNA during cardiac differentiation, we examined the pathways and functions that are associated with their predicated target mRNAs. Therefore, their targets were predicted and extracted for pathway enrichment analysis. GeneGo was employed as the functional annotation system (Fig. 2). Ten highly representative pathways ranked by their respective p-values are shown in Table S5.
On D0 and during early differentiation of D2, D6 (Table S5A, B), 59/ 111 and 43/111 of miRNA targets were involved in TGF, WNT or cytoskeletal remodeling pathway as shown in Fig. 3. Furthermore, targets of up-regulated miRNA on D9 should be preventing the formation of functional cardiomyocytes. Results indicated that Thromboxane A2 pathway signaling is the major pathway targeted by miRNAs in cardiac myocytes (Table S5C Fig. S6).

Targeted genes by different time points of miRNAs
Differentiation of ESCs to cardiac cells is a complex physiological process, accompanied by activation of transcription factors and signal pathways. Given that activation of genes promoting cardiac differentiation requires specific miRNAs, top 10 genes targeted by different time of miRNAs were identified (Table S6).
Tnrc6b and QK showed the most frequency at different time points. Thus, they should be further investigated. As tables shown: Thrc6b, Dnmt3a and Lrrc58 play important roles in maintaining pluripotent stem cell on D0; Tnrc6b, QK and Tmod2 have vital functions during differentiation on D2 and D6; in addition, Fut9, Tnrc6b and Ccnj take parts in forming pure cardiomyocytes on D9.

Networks of genes during cardiac differentiation of mESCs
Dynamically expressed miRNAs regulatory networks were constructed in Fig. 3. On D0 (Fig. 3A), Atxn1 was simultaneously targeted by the highly enriched miR-101a/b-3p, suggesting it may downregulated at the beginning of differentiation. On D2 and D6 (Fig. 3B,  C), CAV2, Col1a1, Col1a2, etc. were targeted by the down-regulated miR-29a-3p and miR-23b-3p, implying that these genes were upregulated during promoting the appearance of the cardiac progenitor cells. On D9 (Fig. 3D), Smad1 was targeted by up-regulated miR-199a-3p, indicating it would preventing the formation of functional cardiomyocytes. To reveal the changes of srRNA expression during differentiation, small RNAs mapping to 5S, 5.8S, 18S, and 28S of M. musculus ribosomal RNA were extracted and characterized separately. Fig. 4A shows the percentage of different kinds of srRNAs in small RNAs. srRNA of 28S accounted for the highest percentage, suggesting that in mESCs, the most srRNAs are 28s. The percentages of total small RNAs mapping to 5S, 5.8S, 18S, and 28S ribosomal RNA were7.3%, 68.1%, 43.0%, and 18.2%, individually.
The length distribution of total redundant srRNAs is shown in Fig. 4B, and two distinct peaks can be observed at nucleotide 19 and 21. Additionally, srRNAs on D2 have the highest counts. These upregulated srRNAs on D2 and D6 suggest that a large number of srRNAs are not likely generated from random degradation [21]. Length distribution of unique srRNAs is displayed in Fig. 4C.
The most abundant 20 srRNAs in four samples are shown in Fig. 4D, and the srRNAs were named according to their position in ribosomal RNA and length [21]. If the position of srRNA is 0, it was mapped before the first nucleotide of the ribosomal RNA. Names and detailed sequences of these srRNAs were shown in Table S7. All of these imply that these srRNAs are up-regulated during differentiation.

Positions of srRNAs in M. musculus ribosomal RNA
Positions of srRNAs in different type of ribosomal RNAs were shown in Fig. 5. In 28s, the higher expressed srRNAs in all samples are the same: 28s-1-21 and 28s-1-19, moreover, they had the similar pattern (Fig. 5A). The normalized reads counts of 28s-1-21 and 28s-1-19 have the same changing pattern during differentiation (Fig. 5B).
To identify genes that are potentially regulated by srRNAs, we performed target analysis for the 20 most abundant srRNAs. The obtained targets were further used in functional ontology enrichment analysis. Ten pathways with lowest p-value are shown (Table S8). Most of these targets categorized into the TGF, WNT and cytoskeletal remodeling pathway (Fig. 3).

Existence and the differential expression of piRNA during differentiation
The PIWI-piRNA pathway has been commonly perceived as germline-specific. Unexpectedly, we found that piRNAs were also expressed during cardiac differentiation. Table S9 lists the top 20 most-abundant piRNAs during differentiation. Most piRNAs were mapped to unique intergenic location, whereas oocyte_piRNA1831, oocyte_piRNA647, oocyte_piRNA3386 have four genome locations, and they are localized in clusters on chromosome 7, 9, 13, 14, respectively. Expression profiling of normalized reads counts were presented. For instance, piR-1078 was mildly expressed on D0, but the expression increased by 9 folds on D2, and reduced gradually. These data, for the first time, suggest that piRNA pathway might be an important regulator for controlling embryonic stem cell maintenance and/or differentiation progression.

Endogenous siRNAs during cardiac differentiation of mESCs
Somatic cells produce siRNAs from exogenous double-stranded RNA (dsRNA) as a defense against viral infection. Detected endogenous siRNAs were displayed in Fig. 6 and Fig. S7. All of these siRNAs duplexes derive from long dsRNA, comprise two~21 nt strands, paired so as to leave 3′, two-nucleotide overhangs. The most abundantly expressed siRNA-0001 is derived from transcriptome SRR1028908.8092784. siRNA sequence, antisense sequence and normalized reads counts are shown in Table S10. All this information implies that endo-siRNAs are not only in C. elegans and flies, but also in mammals, especially in the process of cardiac differentiation from mESCs.

Discussion
For further understanding the regulatory small RNA networks that active in stem cells and cardiac differentiation, we characterized the expression profiling of four stages from mESC to cardiomyocytes.
Using next-generation sequencing, we present a comprehensive analysis of miRNAs in cardiac development. miR-199a-3p was the most abundant miRNA and had the highest expression level on D9. Our results show that NLK is one of its target and NLK was an identified miR-199a-3p targeted gene [6]. NLK could positive effect Wnt signaling pathway, and interact with PPARG leading to transcriptional silencing [12,16]. miR-101a/b-3p was highly enriched on D0, down-regulated on D2, and gradually increased in the following days. Previous studies have shown it could target c-Fos and regulate transforming growth factor-β (TGF) pathway. Our findings showed during early differentiation, targeted genes by down-regulated miRNAs are associated with TGF, WNT and cytoskeletal remodeling pathway. miR-29a-3p was especially down-regulated on D2 and D6. Reports shown that miR-29 family (miR-29a, -29b, and -29c), which negatively regulates Dnmt3a and Dnmt3b, was examined in association with the Wnt/β-catenin signaling pathway [17]. This was also verified in our study. In addition, we also confirmed that the targeted gene of miR-29a-3p was Dnmt3a [2], suggesting our targets prediction is valid.
Targets of specifically expressed miRNAs were assembled in the crucial pathway-TGF, WNT and cytoskeletal remodeling. The TGF and WNT were cross-regulated with epithelial-mesenchymal transition [4]. Wnt-connected transcription factors (Oct-4, Nanog) are primarily responsible for cell proliferation [14]. Moreover, transforming growth factor-β2 has also been shown to enhance differentiation of cardiac myocytes from embryonic stem cells [15]. Thus, this pathway plays a critical role in cardiac differentiation.
In order to identify targeted genes and pathways critical for cardiac differentiation, we focused on the targets of specifically expressed miRNAs. Tnrc6b has been shown to interact with Argonaute proteins, and is required for miRNA-guided gene silencing in HeLa cells [13]. Recent work suggested that Qk contributes to gene expression programs during differentiation [5]. Dnmt3a is known to be expressed in pluripotent stem cells [3].
Our results revealed that, during cardiac differentiation, the expression of srRNAs is up-regulated and inversely related with miRNAs. GeneGo enrichment analysis further revealed that these up-regulated srRNAs targeted genes were associated with TGF, WNT and cytoskeletal remodeling pathway. These findings may provide new insights into the expression regulation of srRNAs in stem cells, and indicate that they could potentially have important yet not been reported biological functions.
Our results provide a global view of RNAome, including miRNAs, srRNAs, piRNAs, endo-siRNAs and tRNAs, as well as their expression pattern during differentiation. It is a useful resource for the researchers to find out the regulatory role of miRNAs and srRNAs and their targeted signaling pathway during the differentiation of ESCs to cardiomyocytes.