Genome-wide profiling and differential expression of microRNA in rat pluripotent stem cells

MicroRNAs (miRNAs) constitute a class of small noncoding RNAs that plays an important role in the post-transcriptional regulation of gene expression. Much evidence has demonstrated that miRNAs are involved in regulating the human and mouse pluripotency. Nevertheless, to our knowledge, miRNAs in the pluripotent stem cells of one of the most commonly used model organisms – the Rattus norvegicus have not been studied. In the present study, we performed deep sequencing of small RNA molecules in the embryonic fibroblasts, embryonic stem cells, and induced pluripotent stem cells of laboratory rats. Bioinformatics analysis revealed 674 known miRNAs and 394 novel miRNA candidates in all of the samples. Expression of known pluripotency-associated miRNAs, such as the miR-290–295 and miR-183-96-182 clusters as well as members of the miR-200 family, was detected in rat pluripotent stem cells. Analysis of the targets of differentially expressed known and novel miRNAs showed their involvement in the regulation of pluripotency and the reprogramming process in rats. Bioinformatics and systems biology approaches identified potential pathways that are regulated by these miRNAs. This study contributes to our understanding of miRNAs in the regulation of pluripotency and cell reprogramming in the laboratory rat.

Small RNA extraction, sequencing and qRT-PCR analysis. Cells were resuspended in RNAlater and used for small RNA extraction with the mirVana ™ miRNA Isolation Kit (Ambion, #AM1560). RNA quality was analysed using an Agilent Bioanalyzer. A 50-ng quantity of small RNA was used to prepare a library, which was then sequenced on the Illumina GAIIX. Raw data were deposited into the SRA database with the accession numbers SRX395433-395442.
Small RNA for cDNA synthesis were extracted from cells and tissues with the mirVana ™ miRNA Isolation Kit. cDNA was synthesized by the NCode ™ VILO ™ miRNA cDNA synthesis Kit (Invitrogen) from 1 µg of total RNA. Real-time PCR was performed using the HS-qPCR SYBR Blue mix (Biolabmix) with a custom miRNA-specific primer and a universal primer from Invitrogen on the LightCycler480 System (Roche). Reaction conditions were 95 °C for 5 min followed by 40 cycles of 95 °C for 15 secs and 60 °C for 30 sec. Reaction specificity was verified by melting curve analysis. Reactions were performed in triplicate. Sequences of the miRNA-specific primers are provided in Supplementary Table S1. Data were normalized to Rnu1a expression.
Mapping and filtration of reads. Raw reads of 11 samples were processed with the FASTX Toolkit v0.0.14 in two steps. First, the data were trimmed with the fastq_quality_trimmer tool, leaving only reads with a tail quality of at least 30. Then, an adapter sequence was clipped with the fastx_clipper program, leaving only 15-nt reads or longer.
Filtered reads were mapped to the R. norvegicus reference genome (version rn5) from the UCSC Genome Browser [29][30][31][32] combined with the ribosomal RNA reference for R. norvegicus from the Silva database 33 . The BWA v0.7.10-r806 software was utilized for gapped/ungapped alignment of single-ended reads. If a read could be mapped to several places in the reference, only one random place was chosen. The portion of reads that mapped solely to chr1-20,X varied from 94.2% to 97.4% in different samples.
To assess the quality of reads, a set of 486 genomic intervals of known premature miRNAs for R. norvegicus was downloaded from miRBase v21 (ftp://mirbase.org/pub/mirbase/CURRENT/genomes/rno.gff3) 34 . The absolute majority (99%) of reads that mapped to the genomic intervals of pre-miRNA had a length of 18-27 nt. To keep the focus of this study on miRNA, we filtered out reads that were shorter than 18 and longer than 27 nt because they likely represented short RNA species, which are different from the miRNA class. The rest of the reads were further filtered from the ones mapped onto parts of the genome that is annotated with rmsk repeats (http:// hgdownload.soe.ucsc.edu/goldenPath/rn5/database/rmsk.txt.gz) 29 of class "RNA" (including rRNA, scRNA, snRNA, tRNA, srpRNA, RNA, and MARNA). This filtration did not influence any reads that were mapped to known miRNAs. A final set of mapped and filtered reads was called informative reads. A summary is presented in Table 1.
Novel miRNA discovery with the custom peaks2mirna method. Wig files were created from the genome maps of the informative reads for two separate groups of samples: PSCs and EFs. The wig files contained RPM (reads per million) expression values for each position on each chromosome. RPM is calculated as the number of reads overlapping a position in the genome divided by the total number of informative reads in a group of samples. Then, wig files were inspected for genome intervals with specific double-peak signals, as explained in Supplementary Figure S1, with a custom Perl script using the following parameters: (i) peak height of 0.9 RPM or higher; (ii) peak length of 18-50 nt; (iii) distance between peaks of 8-30 nt; and (iv) no peaks within 50 nt vicinity. These parameters were selected to maximize the number of detectable, expressed, and known rat miRBase pre-miRNAs (245 of 403, or 60.7%). Then, the following additional filters were applied to the 718 newly identified potential pre-miRNAs to retain only the high-quality candidates: (1) no overlap with the 486 known rat miRBase miRNAs, (2) no overlap with the RNA repeats from the UCSC rn5 rmsk.txt (http://hgdownload.soe.ucsc.edu/ goldenPath/rn5/database/rmsk.txt.gz) 29 , and (3) RNAFold program 35 predicted a hairpin structure. As a result, 73 potential pre-miRNAs remained. Following further analyses of sequence homology, some of them were found to belong to the snoRNA class and had to be analysed separately.
Novel miRNA discovery with the miRanalyzer tool. The miRanalyzer program 36 was applied to each sample and entire dataset of 11 samples. The R. norvegicus genome (assembly version rn4) was used. The default parameters were used except for the "Fast detection of novel microRNAs", which was unmarked. For each sample and separately for entire dataset of 11 samples, potential miRNAs that were named "perfect Dicer" were selected from the predictions. Then, only "perfect Dicers" that did not intersect with genome elements (empty "genome elements" column) were selected. Potential miRNA coordinates in rn4 were converted to rn5 using the UCSC LiftOver tool and then intersected with refSeq genes ("refGene.txt" downloaded from UCSC on 14.04.15). Only the entries that did not overlap the refSeq genes in rn4 and rn5 were kept. If a potential miRNA was found in at least two samples, it was labelled as a "cluster" miRNA. The longest or the most strongly expressed isoform was selected as a representative of the cluster. Potential novel miRNAs that were found by miRanalyzer and were predicted 35 to have a hairpin structure were filtered against known rat miRNAs (ftp://mirbase.org/pub/mirbase/ CURRENT/genomes/rno.gff3) 34 and RNA repeats (http://hgdownload.soe.ucsc.edu/goldenPath/rn5/database/ rmsk.txt.gz) 29 . As a result, 13 potential pre-miRNAs remained.
Search for homologs of novel miRNAs. We took advantage of the UCSC genome browser tracks of rat and other species' genes, including miRNA, that were mapped to the rn6 genome. To find homologs of novel miR-NAs, we intersected their genomic intervals with exons and introns of the ensGene track (last updated 14 Sept, 2015) and the xenoRefGene track (last updated June 19, 2016) (http://hgdownload.soe.ucsc.edu/goldenPath/rn6/ database/) 29 , with consideration of the DNA strand.
Calculation of digital gene expression (DGE) of known and novel miRNA. Coordinates in the rn5 genome of 832 mature miRNAs were downloaded from miRBase v21 (ftp://mirbase.org/pub/mirbase/ CURRENT/genomes/rno.gff3) 34 . The number of All mapped reads overlapping each miRNA genome interval, with consideration of DNA strand, was calculated in each sample. If a read overlapped the border of an interval, it was counted as 1. The sum of overlapping reads served as a raw digital expression level for each miRNA. Alternatively, raw expression was calculated for the informative reads (a subset of All reads that were 18-27 nt in length without RNA-type repeats). The raw expression of known miRNAs calculated by an alternative method was only 1% less than the expression that was calculated based on All mapped reads. To be less stringent, we took the expression level that was based on All reads for downstream analyses. In order to account for the different initial sequencing depths of samples, the RPM normalization procedure was applied to the raw expression values: where i = 1, …, 11 samples; j = 1, …, 832 miRNAs; X i j -raw expression; x i j -RPM expression; and Informative_ reads i -number of informative reads in sample i (Supplementary Table S2). The expression values of novel miR-NAs, based on their genome coordinates, were calculated in a similar way as in the known miRNAs with three exceptions (Supplementary Tables S3 and S4). First, primary sequences of novel pre-miRNAs were mapped to the rn6 genome with the blastn program 38 . If genome intervals of the blastn hits of two different pre-miRNAs overlap, their expression values were summed together. It was done to account for the random locations in the genome for multimapped reads. Second, novel pre-miRNA intervals longer than 50 nt were split into two equal subintervals to mimic possible mature miRNAs in the 5′ and 3′ shoulders. Third, alternative DGE levels were computed using just informative reads to calculate the X i j in formula (1). Both types of DGE based on All and informative reads were used to further classify novel miRNAs to Tier 1 and Tier 2. Tier 1 contained the novel miRNA species that were most likely to be real.

Classification of novel miRNAs into tiers.
By default, all potential pre-miRNAs were located in Tier 1.
For each item, the following quality scores were calculated in order to sort them into two different tiers. First, the ratio of the number of informative reads to the number of All mapped reads was used to assess the portion of 18-27 nt reads. If this ratio was less than 0.5, the novel RNA molecule was moved to Tier 2. If the ratio was between 0.5 and 0.85 and the difference between the counts of All and informative reads was higher than 3, the pre-miRNA was also moved to Tier 2. The second classification parameter was calculated for novel pre-miRNA candidates in which the genomic interval was 50 nt or longer. Such intervals were divided into left and right parts of equal length. The parameter was used to assess the level of expression of the stem loop. It was calculated as the ratio m/(l + r), where l and r are the number of reads overlapping the left and right parts of the interval, respectively (they share reads overlapping the internal border), and m is the number of reads overlapping the border between the left and right parts. If the ratio was higher than 0.05, the potential pre-miRNA was moved to Tier 2.

Principal component analysis (PCA) of miRNA expression levels. RPM expression data of 832
known miRNAs, excluding 158 zero-expression entries, were used. Additionally, 304 novel miRNAs with unique sequences and not homologous to known miRNAs were used independently. Using the R statistical computing program (http://www.R-project.org/) 39 , data were centred around the mean expression values, scaled to unit-variance and then processed using PCA 40 .
Differential expression analysis of miRNAs. ESC and iPSC groups were combined together (PSCs) and compared with EFs. P-values were calculated with the edgeR software 41 for each of the 674 expressed rat miR-NAs from miRBase, assuming an overdispersed Poisson model of raw (not RPM normalized) expression counts. Parameters of overdispersion were fitted by the software, taking into account the common dispersion, trended dispersion and tagwise dispersion of the miRNA expression values. A set of 227 miRNAs satisfied a p-value cut-off of 10 −4 , which corresponded to the 3.0 × 10 −4 false discovery rate (Benjamini-Hochberg). Alternatively, ANOVA was also performed using the 'aov' function in R after replacing zero values with uniform random values from interval (0 to 0.1) and log2 transformation (Supplementary Table S2). However, p-values based on just two samples in one of the groups (EFs) might have little credit. To get over this problem, differentially expressed miRNAs were required to satisfy very strong fold-change conditions without considering the p-value. The ratio between the group averages was required to be at least 3. Additionally, miRNAs upregulated in the PSC group must have a minimum RPM expression in PSC samples that is 2 times (or 10 times in a parallel analysis) higher than its maximum level in the EFs samples. We called this the extremum fold change method. In the same way, miRNAs were defined as upregulated in EFs if their minimum RPM in EFs was 2 times (or 10 times) higher than its maximum level in PSCs, and additionally, they had to have at least three mapped reads in each of the EF samples. Only 219 miRNAs passed this procedure of selection, and they had an edgeR p-value of less than 0.01. Differentially expressed novel miRNAs were defined in the same way.
Clustering according to miRNA expression levels. Expression values of zero were replaced with 0.25, which is below the minimum positive value. After that, log2 transformation was applied. The clustering of miR-NAs and samples based on expression values was performed using the 'pheatmap' function in R using the default parameters.
Search for target genes of differentially expressed miRBase miRNAs. Known rat, mouse and human miRNA-target interactions were downloaded from the MiRTarBase database 42 (http://mirtarbase.mbc. nctu.edu.tw/php/download.php; Catalogued by Species: hsa_MTI.xlsx, mmu_MTI.xls, rno_MTI.xls; accessed on 26/07/2016). Lists of the targets of differentially expressed known miRNAs were obtained taking the miRNA names as the key. They were additionally filtered using gene expression data obtained from sequencing experiments that were performed on the same biological samples 27 . Upregulated targets were kept if the corresponding miRNA was downregulated, and downregulated targets were kept if the corresponding miRNA was upregulated (Supplementary Table S5).
Search for the target genes of novel miRNAs. Mature sequences of novel miRNAs were analysed for targets in the 3′ UTRs of rat refSeq transcripts using the TargetSpy v1.1 43 , miRanda v3.3a 44 and TargetScan v7.0 45 tools. To reduce the false positive rate, only targets that were predicted by all three programs were accepted (Supplementary Table S6).

KEGG and GO enrichment analysis.
Lists of genes found to be targeted by miRNAs were subjected to KEGG and GO enrichment analysis with the DAVID 6.8 Beta online tool 46 using the default parameters. All rat genes were used as the background, and alternatively, only differentially expressed genes 27 were used.

Results
Small RNA sequencing. Small RNA libraries from two EF (RNFM1, RNFF1), two ESC (RES27, dB50), and seven iPSC samples (NF13, NF21, QV8, SU3, one technical, and two biological replicates) were generated and sequenced on the Illumina platform. Reads were filtered as described in the Materials and Methods section. A total of 32.9 million raw reads were generated by the sequencing of eleven libraries followed by trimming and adapter clipping. An absolute majority (31.9 million raw reads) mapped onto the reference genome, including 31.4 million reads that mapped on the chr1-20, X (All reads). After filtration against RNA-type repeats and after keeping only 18-27 nt sequences, a set of 8,860,854 informative reads was obtained. The exact amount of informative reads for each library is presented in Table 1. The length distribution of the informative reads was similar among different samples. Examples are given in Fig. 1a and b. Most of the informative reads were between 21-25 nt. To identify known miRNAs, coordinates of All reads that were mapped to the genome were compared with known miRNA coordinates from miRBase v21 34 . The portion of known miRNA reads (Table 1) varied from 30.84% (dB50) to 81.78% (RNFM1). In total, 674 known miRNAs were identified among all samples. Novel miRNA discovery. To discover novel miRNA species, the informative reads that were mapped to the rn5 genome were processed by three different methods: the custom peaks2mirna, miRanalyzer 36 and miR-Deep2 37 . Pre-miRNAs discovered by more than one method were collapsed by taking the maximum of their genomic intervals. Likely false positives, in which the primary sequence was repeated in the genome more than 15 times, were excluded. The number of repeats in the genome was calculated as the number of blastn 38 hits, and the threshold was fitted to avoid false negatives. Finally, 388 potential novel pre-miRNAs including genomic duplicates were identified. This number included 39 pre-miRNAs that were homologous to rat or non-rat miR-NAs and 40 molecules that were homologous to small nucleolar RNA (snoRNA). Homology analysis was done by overlapping genome intervals with the xenoRefGene and ensGene UCSC genome tracks 29 . This analysis allows nucleotide gaps, deletions and mismatches between homologs, including variance in miRNA seed sequence. The rest of the 309 candidates were split into Tier 1 and Tier 2 lists according to their quality scores as described in the Materials and Methods section and each had 248 and 61 items, respectively. Tier 1 and novel species that are homologous to known miRNAs are presented in Supplementary

Differential expression of miRNAs in rat pluripotent stem cells and embryonic fibroblasts.
To reveal differentially expressed miRNAs, reads per million (RPM) gene expression was calculated for the 832 known and the 557 potential novel miRNAs using All mapped reads as described in the Materials and Methods section. Among the known miRNAs, 77 were found to be upregulated and 142 were found to be downregulated in PSCs compared to EFs, with at least a double extremum fold change (Supplementary Table S7). For novel Tier 1 miRNAs, the corresponding up-and downregulated genes were 10 and 13, respectively (Table 2). In Tier 2, there were 9 and 2, respectively (Supplementary Table S4). The criteria for miRNA up-or downregulation were described in the Materials and Methods section. Most of the miRNAs that were upregulated in PSCs have expression levels that are below the median in adult rat tissues, according to RT-PCR expression profiling with the TaqMan ® Rodent miRNA card A v2.0 (GEO accession GSE66195, https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE66195). On the other hand, most of the downregulated miRNAs have high expression levels (above median) in somatic tissues in the same experiment (data not shown). PCA was performed based on the RPM values of 219 differentially expressed known miRNAs. Figure 1c shows the PCA projection of eleven cell lines onto the first two principal components (PC). The first PC perfectly separated EF samples from the others. The second PC divided samples by gender into two groups: NF13, NF21 and RES27 were female and SU3, QV8 and dB50 were male. PC1 explained 81.1% while PC2 explained 6.6% of the total variance. PC1 also perfectly separated PSCs from EFs if considering all (not just differentially expressed) known or novel miRNAs (data not shown). We then performed sample clustering based on differentially expressed miRBase miRNAs (extremum fold change of at least 10) and on all expressed miRNAs. Figure 1d and e show heat maps of the clusters of miRNAs that are overexpressed in PSCs and EFs, respectively. There was no significant difference between the results of the hierarchical clustering of samples performed on all miRNAs and on only differentially expressed miRNAs (data not shown). Both methods, PCA and hierarchical clustering, revealed that there are principal differences in miRNA expression between EFs and PSCs but no major differences between ESCs and iPSCs.
A number of miRNAs that were found to be upregulated in rat PSCs, such as the miR-290-295 cluster and the miR-200 family miRNAs, are involved in pluripotency induction and maintenance in mice 21,47 . Several miRNAs, such as members of the miR-183-96-182 cluster, miR-741, miR-881, miR-878, miR-743a, miR-743b, miR-465, miR-871, miR-880, are known to be expressed in mouse ESCs and were also expressed in rat PSCs; however, their function in pluripotency induction and regulation has not been directly shown 48-50 . GO and pathway analyses of targets of the differentially expressed miRNAs. Analysis of known miRNA targets. To reveal the function of differentially expressed miRNAs in rat PSCs and EFs, we performed  an analysis of the published targets of these miRNAs in humans, mice, rats, and other species. We took known miRNAs with an extremum fold change of more than 10 between EFs and PSCs. Among them, 58 miRNAs were upregulated in the PSCs and 69 were upregulated in the EFs (Supplementary Table S7). Known miRNA targets were found using information from the MiRTarBase database 42 . Additionally, the targets were required to be downregulated for upregulated miRNAs and upregulated for downregulated miRNAs, which is consistent with the expected miRNA function. By this method, we found 388 published target genes for miRNAs that were upregulated in PSCs and 69 published targets for downregulated miRNAs (Supplementary Table S5). All of the upand downregulated target genes were analysed with the DAVID 46 web application tool to identify the associated KEGG pathways and GO terms. KEGG pathways enriched with target genes, using all rat genes as a background, are presented in Fig. 2a and in Supplementary Table S8. Most of the target genes were associated with cancer pathways, stem cell pluripotency regulatory pathways, and the Hippo signalling pathway. Of note is that the results persist if only differentially expressed genes were used as a background (data not shown). Among the enriched GO terms, stem cell differentiation, response to retinoic acid, somatic stem cell population maintenance, regulation of apoptotic process were present (Fig. 2b). The obtained results suggest that these differentially expressed miRNAs play a substantial role in the regulation of pluripotency and cell reprogramming in rats.
Analysis of novel miRNA putative targets. Targets of 304 novel miRNAs with unique sequences that are not homologous to known miRNAs were predicted by TargetSpy v1.1 43 , miRanda v3.3a 44 , and TargetScan v7.0 45 as described in the Materials and Methods section. In total 8,437 transcript targets were identified for novel miRNAs (Supplementary Table S6). GO enrichment and KEGG signalling pathway analyses were performed using the DAVID tool 46 on the targets of the top 10 highly expressed novel miRNAs in the PSC group and the top 10 highly expressed novel miRNAs in the EF group, including all rat genes as a background. The results are presented in Fig. 3 and Supplementary Table S9. According to the GO annotations, the predicted target genes of novel miRNAs expressed in PSCs are related to processes of embryonic development, cell migration, signal transduction and gene expression regulation. According to KEGG annotation, putative targets are involved in endocytosis, pathways in cancer and the Hippo, FoxO, Wnt, and TGF-beta signalling pathways (Fig. 3a,b). The two latter networks are known to be involved in cell reprogramming and the regulation of pluripotency in rats 12,47 .

Validation of novel miRNAs by qRT-PCR.
To validate some of the novel miRNAs, we performed qRT-PCR analysis using the NCode ™ system in the EF and PSC lines. Four novel miRNAs with the following coordinates: chrX:−:151288045-151288101 (rno-miR-novel-8); chr7:+:70463555-70463594 (rno-homolog-miR-26); chr3:+:16697111-16697143 (rno-homolog-miR-199); and chr18:−:69422790-69422857 (rno-miR-sno-57) were selected for the validation. rno-miR-novel-8, rno-homolog-miR-26, and rno-homolog-miR-199 miRNAs were selected from Tier 1, and rno-miR-sno-57 miRNA was selected from Tier 2. In addition, we analysed the expression of miR-741-3p and miR-743a-3p and found that, in accordance with sequencing data, they were highly expressed in rat PSCs. A pluripotency-associated miRNA, miR-295-3p, was used as a positive control. Expression of miR-295-3p was observed in rat PSCs and absent in rat EFs, which was consistent with its expression pattern in the mouse 51 (Fig. 4a). qRT-PCR analysis demonstrated that the rno-miR-novel-8 and rno-miR-sno-57 miRNAs as well as the known miRNAs miR-741-3p and miR-743a-3p were expressed in rat PSCs and absent in rat EFs (Fig. 4a). rno-homolog-miR-26 and rno-homolog-miR-199 miRNAs were expressed in EFs, ESCs and iPSCs, which is consistent with the data obtained from sequencing. We also analysed the expression of the rno-miR-novel-8 miRNA, rno-miR-sno-57 miRNA, miR-741-3p, and miR-743a-3p, which are upregulated in PSCs during differentiation. We induced spontaneous differentiation of rat ESCs and iPSCs in a monolayer and analysed the expression of the miRNAs at various time points. We found that rat miR-295-3p, miR-741-3p and miR-743a-3p expression was significantly decreased during ESC and iPSC differentiation (Fig. 4b). For novel miRNAs, the expression of rno-miR-novel-8 miRNA was absent at the first time point (day 4) and subsequent time points of differentiation, whereas the expression of rno-miR-sno-57 miRNA increased at the first stages of differentiation and then decreased to levels lower than that in PSCs (Fig. 4b). We also assessed the expression of rno-miR-novel-8 miRNA, rno-miR-sno-57 miRNA, and the known miRNAs in different rat adult organs, including the heart, lungs, liver, kidneys, ovaries, testes, eyes and spleen. We found that miR-741-3p, miR-743a-3p, and rno-miR-novel-8 miRNA were expressed in the testis at levels comparable to that in PSCs. miR-741-3p was also expressed in the lungs, ovaries and spleen but at lower levels (Fig. 4c).

Discussion
R. norvegicus is one of the most promising models for disease modelling and translational biomedicine. Pluripotent stem cells are currently widely used to generate various disease models 52 . An understanding of pluripotency and self-renewal regulation is required to improve the approaches to PSCs generation and cultivation, which will help realize the full potential of their application. Although much progress has been made in the understanding of the pluripotency regulatory network, it has not been fully investigated yet, especially in rats. miRNAs are known to play substantial roles in the regulation of pluripotency and self-renewal 14 . To study miRNA expression in rat PSCs and to find novel rat miRNAs, we performed genome-wide sequencing of small RNA libraries from rat ESCs, iPSCs, and EFs. Analysis of the sequencing data revealed the expression of 674 known miRNAs across all samples. Differential expression analysis of known miRNAs demonstrated that 77 miRNAs were upregulated and 142 miRNAs were downregulated in rat PSCs compared to EFs. Members of the miR-290-295 cluster (miR-290, miR-291a, miR-292, miR-291b, miR-293, miR-294 and miR-295-1) were the most abundant among known miRNAs that were upregulated in rat PSCs. It is known that mouse miR-290-295 and its human homologues in the miR-372 clusters are ESC specific and regulate cell cycle progression 20,21 . Among the miRNAs upregulated in rat PSCs, we also identified miR-205 and members of the miR-200 family (miR-200a, miR-200b, miR-200c, miR-141 and miR-429), which promote mesenchymal to epithelial transition (MET) in mouse cells, a key step in fibroblast reprogramming 47 . Members of the miR-200 family are known to suppress the Zeb2 transcription factor, and its inhibition promotes MET 53 . Another miRNA cluster, miR-183-96-182, is known to be expressed in mouse PSCs 49 and is also expressed in rat PSCs. The microRNA-183-96-182 cluster inhibits epithelial to mesenchymal transition (EMT) in cancer cells, and thus, the miRNA cluster may facilitate the reprogramming process by promoting MET [54][55][56][57] . Members of the miR-106-363 and miR-106b-25 clusters, which were shown to be upregulated during fibroblast reprogramming in mice 58 , were also upregulated in rat PSCs. Most of the miRNAs upregulated in rat PSCs are also expressed in mouse ESCs or iPSCs 48 . In addition, we found several rat miRNAs (including miR-3580, miR-3551, and miR-6325) whose expression had not been established in PSCs before.
A cluster of miRNAs, comprised of miR-463, miR-465, miR-471, miR-741, miR-743a, miR-743b, miR-871, miR-878, miR-880, miR-881, and miR-883, was located on the X chromosome and was expressed at a high level in rat PSCs according to the sequencing data. These miRNAs are known to be expressed in mouse ESCs and testes 48,50,59 . These miRNAs are testis specific in rats, according to the miRNA sequencing data of 21 male rat organs 26 . Using qPCR, we have shown that members of the miR-741-3p and miR-743a-3p clusters were expressed at high levels in rat ESCs and iPSCs and were downregulated during differentiation. These miRNAs are also expressed in some adult rat tissues, according to our data and RT-PCR expression profiling with the TaqMan ® Rodent miRNA card A v2.0 (GEO accession GSE66195, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66195).
The number of known rat miRNAs is much lower compared to that in humans and mice. In this study, 394 novel mature miRNAs were identified using three methods -mirDeep2, miRanalyzer, and custom peaks2mirna. Most of the novel miRNAs were identified by the mirDeep2 method. In rat PSCs and EFs, 23 novel miRNAs were differentially expressed. Among them, 10 novel miRNAs were upregulated and 13 were downregulated in PSCs. This list includes homologues of human miR-302a and miR-302b. We found that these miRNAs were upregulated in rat PSCs. We also found a homologue of human miR-367. The human and mouse miR-302-367 cluster has been demonstrated to play an important role in the regulation of pluripotency. Moreover, overexpression of the miRNA cluster was sufficient to reprogramme fibroblasts to iPSCs 21,23,60 . However, expression levels of miR-302a and miR-302b were low in rat PSCs compared to the ESC specific miR-290-295 cluster. This is in agreement with the fact that high expression levels of miR-302 cluster are typical of mouse epiblast stem cells 48 . We also found a homologue of human mir-372 located inside the rat miR-290-295 cluster. miR-372 in humans is known to be connected to fibroblasts reprogramming, particularly by promoting MET 60 .
Four novel miRNAs were validated using qRT-PCR analysis. Three miRNAs were selected from Tier 1, and one was selected from Tier 2. Expression of the two validated novel miRNAs was observed in rat PSCs and was absent in EFs. Expression of rno-miR-novel-8 miRNA decreased during differentiation and was not observed in the analysed adult organs, except in the testes. The pattern is similar to ESC-specific miRNAs, such as members of the miR-290-295 cluster 51 . In addition, we found that the rno-miR-novel-8 miRNA was located less than 10 kb from miR-871, mir-3580 and miR-465, which were upregulated in rat PSCs. This may be a cluster of miRNAs that is involved in the regulation of pluripotency, and further investigations are required to shed light on its function. miRNAs implement their functions through repressing target gene expression. In human and mouse PSCs, clusters of pluripotency-associated miRNAs control cell proliferation and prevent differentiation 58,61 . To better understand the functions of known rat miRNAs that are differentially expressed between PSCs and EFs, we analysed their target genes that have been validated in rats, humans or mice. An additional criterion was that the change in the PSC and EF target gene expression levels had to be opposite of that for miRNAs. This criterion helped to improve the analysis of miRNA functions in the rat. In total, 388 downregulated targets were found for the upregulated miRNAs in PSCs, and 69 upregulated targets were observed for miRNAs that were downregulated in PSCs. Then, we performed KEGG and GO enrichment analyses of the target genes. The most common GO terms were related to transcriptional regulation, which indicated that a substantial part of the target genes functioned as transcriptional factors. Many GO terms, such as stem cell differentiation and response to retinoic acid, were related to cell differentiation, which was consistent with the function of miRNAs in suppressing PSCs differentiation 61, 62 . Self-renewal and differentiation processes are balanced in PSCs, and miRNAs are involved in the regulation of this process 58,61 . We also found a stem cell population maintenance category that was relevant for targets of miRNAs that were differentially expressed in PSCs. Regulation of the apoptotic process was also an important category for pluripotent cells, because inhibition of apoptosis promotes pluripotency 63 and was demonstrated for the ESC-specific miR-290-295 cluster 64 . We found enrichment of differentially expressed miR-NAs and their targets in the KEGG signalling pathway, the stem cell pluripotency regulatory pathway, pathways in cancer, and the Hippo signalling pathway. Pathways in cancer are relevant for miRNAs that are upregulated in PSCs. It is known that miRNAs that are expressed in PSCs are frequently related to different types of cancer. The human miR-372-373 cluster promotes tumorigenesis, for example 65 . The Hippo signalling pathway represents a barrier for reprogramming, and its inhibition by LATS2 knockdown increases the efficiency of iPSC generation 66 . Analysis of the putative target genes of novel rat miRNAs has shown that the miRNAs may take part in the regulation of development and differentiation. In addition, KEGG signalling pathway analysis suggests the involvement of miRNAs in the regulation of the Wnt signalling pathway, which regulates rat PSC maintenance, and the TGF-beta pathway, which participates in the regulation of cell reprogramming 12,47 . Based on the results, we suggest that the rat miRNAs identified in this study are involved in the regulation of pluripotency and reprogramming in rats.
In conclusion, we first compared the miRNA expression in rat embryonic stem cells and induced pluripotent stem cells with rat embryonic fibroblasts used for reprogramming. The miRNA expression signature in rat PSCs was similar to that in mouse PSCs, and most mouse pluripotency-associated miRNAs were also found in rats. A number of novel rat miRNAs were identified. Overall, our results shed light on the role of miRNAs in the pluripotency regulatory network of the laboratory rat.