Long noncoding RNAs sustain high expression levels of exogenous octamer-binding protein 4 by sponging regulatory microRNAs during cellular reprogramming

Long noncoding RNAs (lncRNAs) modulate gene expression as competing endogenous RNAs (ceRNAs) that sponge regulatory microRNAs (miRNAs). During cellular reprogramming, genes associated with pluripotency establishment need to be up-regulated, and developmental genes need to be silenced. However, how ceRNAs control cellular reprogramming still awaits full elucidation. Here, we used doxycycline-inducible expression of the four transcription factors octamer-binding protein 4 (OCT4), SRY-box 2 (SOX2), Krüppel-like factor 4 (KLF4), and proto-oncogene c-Myc (c-Myc) to generate induced pluripotent stem cells (iPSCs) from mouse embryonic fibroblasts (MEFs). Using RNA-Seq and bioinformatics approaches, we found that the expression levels of miRNAs from MEFs remain high from day 0 to 6 after the doxycycline induction. Many genes targeted by these miRNAs were up-regulated, and long intergenic noncoding RNAs (lincRNAs) and circular RNAs (circRNAs), which have complementary binding sites to these miRNAs, were highly expressed, indicating lincRNAs and circRNAs may function as ceRNAs. Intriguingly, knockdown of the linc/circRNAs that sponge the miRNAs, which target OCT4 down-regulated exogenous OCT4, decreased reprogramming efficiency, and resulted in low-grade iPSCs. Our results suggest that the ceRNA network plays an important role in cellular reprogramming.

The microRNAs (miRNAs) 2 (ϳ22 nucleotide noncoding RNAs) are an abundant class of small noncoding RNAs, which guide the RNA-induced silencing complex to miRNA response elements (MREs) to its targeting transcripts, and can induce mRNA degradation or translational repression (1)(2)(3). The efficiency of miRNAs is known related to the abundance of MREs present in the transcripts (4)(5)(6). Previous studies have shown that one gene can have multiple MREs for distinct miRNAs, and vice versa, miRNA can target on multiple distinct transcripts (7,8). A highly expressed MRE-containing transcript can compete for multiple shared miRNAs and lead to observable changes in miRNA activity and thus regulate the expression of one or multiple targeted transcripts. This "sponge" mechanism has been proposed for competing endogenous RNAs (ceRNAs) to regulate other RNA transcripts by competing for shared miRNAs (9,10).
Such a sponge mechanism involves coding and noncoding transcripts, including long intergenic noncoding RNAs (lin-cRNAs) (10 -13). For example, linc-MD1, a muscle-specific ceRNA, regulates muscle differentiation (13). Recently, another naturally occurring family of noncoding RNAs, the circular RNA (circRNA), from back-spliced exons, has also been reported to act as miRNA sponges to regulate gene expression (14 -17), such as the miR-7 inhibitor in the central nervous system (15). Although thousands of lincRNAs have been annotated in eukaryotic genomes and targeted by miRNAs with both computational and experimental evidences (18,19), so far, only a small number of lincRNAs have been identified to serve as ceRNAs. The circRNAs are also found abundant in eukaryotic transcriptome (20,21), whether circRNAs can function as miRNA sponges remains completely elusive. Moreover, the biological relevance of ceRNAs has recently been challenged, because the relatively low expression level of most lncRNAs might limit their ability to effectively modulate, in a miRNAdependent manner, mRNA abundance. For example, the level of one transcript, Aldoa, required to significantly alter the level of one highly abundant miRNA, miR-122, and its targets in adult hepatocytes was found to exceed the changes observed in vivo, even under extreme physiological or disease conditions (5). The discrepancies in the conclusions of the different attempts suggest that substantially more genetic and genomic evidence from diverse cell types will be required to resolve this issue and to establish the general prevalence and physiological relevance of ceRNAs.
Considering that many established ceRNAs neither share an unusual high number of predicted MREs with their mRNA targets nor are especially abundant (10,22), the ability to modulate mRNA abundance was thought to be limited (4,6). However, accumulating evidences suggest that ceRNAs might play an important role in regulating key transcription factors during the cell-fate decision process, supported by the observation that the expression of ceR-NAs is tightly regulated with spatial and temporal expression patterns. For example, linc-RoR competes for miR-145 binding with key self-renewal transcription factor transcripts, including Nanog, Oct4, and Sox2, and is expressed in undifferentiated embryonic stem cells (ESCs) (11). The specific expression profiles of these ceRNAs might specify the cells in which their activity exerts the greatest effect. Cell-fate decision often involves switch-like responses in the expression levels of key regulatory genes that result in coordinated changes in transcription profiles, suggesting that lowly abundant, but specifically expressed, ceRNAs might function efficiently on regulation of the transition from differentiated to pluripotent cell states.
Cell fate conversion can be induced by over-expressing sets of regulatory transcription factors (23)(24)(25). The over-expression of Oct4, Sox2, Klf4, and c-Myc (OSKM) can reprogram mouse embryonic fibroblasts (MEFs) to induced pluripotent stem cells (iPSCs) (26,27). In this study we performed ribominus RNA-seq and small RNA-seq from cells collected at the early reprogramming stages and fully reprogrammed iPSCs, which can generate "all-iPSCs" pups ( Fig. 1A), to profile expression patterns of miRNAs, lincRNAs, circRNAs, and mRNAs. We also took advantage of publicly available RNA-seq and small RNA-seq data from Mbd3 flox/Ϫ (a member of nucleosome remodeling and deacetylase co-repressor complex) reprogramming system, which results in near 100% efficiency of iPSCs reprogramming within 7 days by OSKM induction (28), to confirm our analysis. Our studies ought to determine the possible function of linc/circRNAs as miRNAs sponge during cellular reprogramming.

Profile of miRNAs from MEFs during reprogramming
MEFs, from mice carrying a single Dox-inducible OSKM polycistronic 4F2A cassette (tetO-OSKM), were induced by reprogramming by adding doxycycline to the culture medium. To demonstrate the miRNA-mediated cross-talk during reprogramming, we first profiled small RNAs in MEFs, the reprogramming MEFs (rMEFs) at day 3 and day 6 post-OSKM induction and the fully reprogrammed iPSCs. To confirm our data, profiling of small RNA (GEO accession GSE102518) from the Mbd3 flox/Ϫ reprogramming system were also analyzed. 451 and 474 miRNAs were co-expressed among MEFs, rMEFs, and iPSCs in our and Mbd3 flox/Ϫ reprogramming systems, respectively, indicating miRNAs from MEFs are sustainable during reprogramming (Fig. 1B). To figure out the extent of miRNAmediated interactions between their mRNA targets and lncRNAs, we selected 311 and 279 miRNAs highly expressed in our and Mbd3 flox/Ϫ MEFs (top 50% abundance and RPM Ͼ20), and 244 miRNAs were overlapped (Fig. 1C). The expression levels of the MEF-high miRNAs were not significantly changed over the course of reprogramming, and significantly decreased in reprogrammed cells (at days 7 and 8 in the Mbd3 flox/Ϫ reprogramming system) and iPSCs (Fig. 1D, Table S1), and the miR-NAs were chosen for further analysis.

Activity of miRNAs may be inhibited during reprogramming
During reprogramming, genes associated with pluripotent establishment need to be activated or up-regulated, whereas developmental genes need to be silenced (28). To determine the role of miRNAs in regulating the transcriptomic changes, we profiled expression of mRNAs on MEFs, rMEFs, and iPSCs as well, and then analyzed the effect of the MEF-high miRNAs on up-and down-regulated genes during reprogramming. We found the miRNAs were not only targeted on the down-regulated genes, but also the up-regulated genes, including key pluripotent transcriptional factors, such as Oct4, Klf4, Nanog, Tbx3, Sall4, and Esrrb. Notably, the number of up-regulated genes targeted by the miRNAs was more than that of the downregulated genes ( Fig. 2A; Table S2). The expression levels of miRNAs targeting on up-and down-regulated genes had no significant difference (p Ͼ 0.05), and both were maintained at high levels during reprogramming and decreased in iPSCs (Fig.  2B). Then, we checked the densities of MREs in up-and downregulated mRNA targets. We found that the MRE density in down-regulated mRNA targets (mean of 0.430 MREs/kb) is not significantly different with up-regulated mRNA targets (mean of 0.487 MREs/kb; p Ͼ 0.05) (Fig. 2C). Considering that highaffinity MREs would be more effective to bind miRNAs than low-affinity MREs, resulting in down-regulation of gene expression more efficiently, we analyzed the ratio of 8-mer, 7-mer, and 6-mer MREs in the targets, and found no significant difference between the up-and down-regulated targets (Fig.  2D). Our results show that neither the density nor the affinity of MREs is different between the miRNA-targeted up-and downregulated genes. During reprogramming, developmental genes should be silenced, so the highly-expressed miRNAs targeting the down-regulated genes may have little effect on reprogramming. However, if expression of the up-regulated genes is influenced by the miRNAs, the reprogramming will fail. Thus, we hypothesize that additional mechanisms may exist to inhibit the activity of the MEF-high miRNAs during reprogramming.

LincRNAs and circRNAs sponge miRNAs during reprogramming
binding sites with the miRNAs, and found seven linc/circRNAs were highly expressed in rMEFs at days 3 and/or 6, compared to MEFs (Fig. 3E). The results suggest the general function of linc/ circRNAs as sponges on inhibition of activity of miRNAs during reprogramming.

Activity of miRNAs targeting on Oct4 is counteracted during iPSCs reprogramming
To further demonstrate the ceRNA mechanism during reprogramming, we tested the hypothesis in the context of Oct4, a key pluripotent transcriptional factor that determines the successful reprogramming and is not replaceable in induction of iPSCs (31). We predicted seven MEF-high miRNAs targeting on the CDS region of Oct4 (Oct4-miRNAs) by RNA22 (32) (Fig. 4A). However, definitely, Oct4 was up-regulated ( Fig.  4B). During iPSCs reprogramming, endogenous Oct4 activation generally takes place at 7-12 days (33), and, in the study, we found that the expression of endogenous Oct4 was initiated from day 7 post-OSKM induction (Fig. 4C), so the ceRNA mechanism we checked mainly involves the regulation of exogenous Oct4 expression. During reprogramming, the expression of the seven Oct4-miRNAs were sustained in high levels, and significantly decreased in reprogrammed cells (at days 7 and 8 in Mbd3 flox/Ϫ reprogramming system) and iPSCs (Fig. 4D). This result suggests that the expression of Oct4 is not tightly regulated by the miRNAs or a mechanism that can counteract the function of the miRNAs. To confirm if Oct4-miRNAs can down-regulate the expression of exogenous Oct4 during the process of reprogramming, the increasing concentrations of the mimics of the miRNAs were transfected into MEFs and the expression of Oct4 mRNA was monitored by qRT-PCR at day 3 post-OSKM induction. We found that mmu-miR-130b-5p, mmu-miR-421-3p, mmu-miR-497a-5p, and mmu-miR-532--3p could efficiently down-regulate Oct4 expression at high concentrations, and the dose-dependent manner revealed that activity of the miRNAs might be limited (Fig. 4E). To ensure that the miRNAs can definitely function on Oct4 mRNA, we inserted the native and mutant MREs of the four Oct4-miRNAs into the 3Ј end of a standard luciferase reporter after an internal ribosome entry site (Fig. 4F) and transfected into MEFs, respectively. The luciferase assay was detected at day 3 post-OSKM induction and showed that the luciferase activity expressed by the reporters with native MREs was significantly decreased, and by the reporters with mutant MREs had no significant changes, suggesting the miRNAs can indeed repress the expression of Oct4 by targeting on the CDS (Fig. 4G). These results suggest that the Oct4-miRNAs can down-regulate exogenous Oct4 expression, however, their activities are counteracted during reprogramming.

Linc/circRNAs can sustain high expression of exogenous Oct4 by sponging the miRNAs targeting on Oct4
We considered if the linc/circRNAs, which have complementary binding sites with the Oct4-miRNAs could counteract

LincRNAs and circRNAs sponge miRNAs during reprogramming
the activity during reprogramming. We indeed found that some linc/circRNAs, which have complementary binding sites with the four Oct4-miRNAs were highly expressed during reprogramming, and decreased in reprogrammed cells and iPSCs (Fig. 5A). We postulate that if we knockdown the expression of the linc/circRNAs, exogenous Oct4 could be down-regulated for the increased activity of the miRNAs. We selected the efficient locked nucleic acid (LNA)-siRNA for each lincRNA and circRNA (LNA-siRNAs targeting the back-splice sequences of circRNAs were used) (Fig. S1). As expected, Oct4 expression was significantly down-regulated at least 20% when some linc/circRNAs were knocked down (Fig. 5B). We supposed that these linc/circRNAs might benefit Oct4 expression through sponging the Oct4-miRNAs. To test the hypothesis, we conducted RNA immunoprecipitation (RIP) using day 3 rMEFs from MEFs transfected with the biotinlabeled miRNAs at a final concentration of 20 nM, and observed that ENSMUSG00000092341, mmu_circ_00006895, and mmu_circ_00000319 were specifically enriched by qRT-PCR analysis normalized to captured Oct4 mRNA (up to 20-fold enrichment compared with Gapdh mRNA), suggesting the three linc/circRNAs are able to interact with the Oct4-miRNAs (Fig.  5C). Further analysis showed multiple complementary binding sites with mmu-miR-130b-5p, mmu-miR-421-3p, mmu-miR-497a-5p, and mmu-miR-532-3p on the sequences of the linc/cir-cRNAs (Fig. S2). LncRNAs in the cytoplasm may function as miRNA sponge. Thus, we analyzed the subcellular localization of the linc/circRNAs. Gapdh and Xist were checked as control. The results showed that the linc/circRNAs were mainly located in the cytoplasm (Fig. 5D). These results show that the linc/circRNAs function as Oct4-miRNAs sponges to sustain high expression of exogenous Oct4 during reprogramming.

The linc/circRNAs can improve reprogramming efficiency and the quality of iPSCs
The efficiency to derive alkaline phosphatase (AP)-positive clones was about 1.31% using the secondary reprogramming system (Fig. 6A), and the iPSCs could give rise to chimeric pups  and contributed to germline transmission (Table 1). Using the reprogramming system, we investigated whether the linc/cir-cRNAs have any functional effects on reprogramming. To this end, we transfected each of the LNA-siRNAs targeting the linc/ circRNAs in MEFs, respectively, and we did not find any significant changes about the formation of AP-positive colonies at day 17 post-OSKM induction. However, when the mixture of the LNA-siRNAs (si-mix) were transfected, we found that colony formation was significantly decreased (Fig. 6A) and only two cell lines were obtained. When additional Oct4 was added back (si-mixϩOct4-addback) by infecting MEFs using lentiviral Oct4, the number of AP-positive colonies was greatly enhanced (Fig. 6A). To test if the overexpression of these lncRNAs can enhance the reprogramming, we cloned a part of the sequence (3000 bp) of ENSMUSG00000092341 (the sequence includes 3, 2, 2, and 1 complementary binding sites of mmu-miR-532-3p, mmu-miR-497a-5p, mmu-miR-421-3p, and mmu-miR-130b-5p, respectively) and overexpressed it in MEFs, and then induced the reprogramming by adding doxycycline. However, after day 17 post-OSKM induction, we did not observe any significant changes in the formation of APpositive colonies (Fig. S3). To further confirm that low reprogramming efficiency is related to the decreased expression of exogenous Oct4 induced by knockdown of the linc/circRNAs, we tested the expression of the other three reprogramming factors, Sox2, Klf4, and c-Myc, after knockdown of the three linc/ circRNAs, and found their expressions were not significantly affected (Fig. S4). Low levels of exogenous Oct4 leads to gener-

LincRNAs and circRNAs sponge miRNAs during reprogramming
ation of iPSCs with aberrant DNA methylation of the Dlk1-Dio3 locus, low expression of miRNAs encoded by the locus, and the low capacity in chimeric mice (33). Thus, we checked the methylation status of the Dlk1-Dio3 locus. Control iPSCs showed normal, 50 -60% DNA methylation, whereas the linc/ circRNA-knockdowned cells, especially the si-mix cells, showed DNA hypermethylation, and Oct4-addback could rescue the methylation defect (Fig. 6B). Consistent with the hypermethylation pattern, the miRNAs in the locus showed significantly lower expression in linc/circRNA-knockdowned cells compared to control iPSCs. Also, the expressions could be rescued by Oct4-addback (Fig. 6C). In addition, we tested the in vivo developmental potency of the iPSCs by injection into blastocysts. In contrast, only one si-mix cell line could produce chimeric mice with the low extent of chimaerism and fail to support germline transmission (Table 1), whereas, all three Oct4-addback cell lines produced chimeric pups and one cell line could contribute to germline transmission (Fig. 6D, Table  1). The results indicate that the linc/circRNAs can improve reprogramming efficiency and the quality of iPSCs by sponging miRNAs targeting on Oct4.

Discussion
Herein we investigated the prevalence and properties of lncRNAs as ceRNA during reprogramming by experimental and bioinformatic approaches. We focused on the analysis of linc/circRNAs with the proposed role in the regulation of cell states from differentiation to pluripotency. Our data reveal the linc/circRNAs can sponge the miRNAs targeting on Oct4 mRNA and regulate expression of exogenous Oct4, and further
In the context of this critical cellular transition when repertoires of miRNAs are limited, the impact of linc/circRNAs on transcriptional programs is likely to be greater than in differentiated cells. Environmental or cellular stress, for example, upon starvation or infection, may also offer similar opportunities for strong effects of ceRNA. We integrated the data of ribo-minus RNA-seq and small RNA-seq from cells collected during reprogramming and fully reprogrammed iPSCs. During reprogramming, OSKM induction causes drastic changes of transcriptome (23,27,34). It leads to inactivation of genes responsible for the identity of the parental cells and the activation of genes that are crucial for the establishment of the lineage of interest (23). Similarly, in our study, we observed initial silencing of fibroblast genes and their transcriptomic activity was gradually replaced with switching-on genes associated with pluripotency. However, we found the pool of miRNAs from MEFs did not changed remarkably during reprogramming, indicating reprogramming of the miRNA expression pattern is a latter event. To ensure the successful reprogramming process, the activity of miRNAs targeting on genes that need to be activated or up-regulated must be reduced. Indeed, we identified many linc/circRNAs, which share miRNA-binding sites with mRNAs of pluripotent transcription factors, and were highly expressed, suggesting that ceRNAs could be beneficial for reprogramming via a miRNAmediated mechanism.
LncRNAs are proposed to act as sponges for miRNAs to modulate gene expression (4,35). In mouse ESCs, many lncRNAs have been predicted to compete for miRNAs with mRNAs (4). However, it has not yet been experimentally validated as ceRNAs, especially circRNAs. A previous study has shown that linc-ROR regulates Oct4, Nanog, and Sox2 to sustain self-renewal of human ESCs by sponging miRNAs (11). In the study, we identified that linc/circRNAs could serve as ceR-NAs during reprogramming. We evidenced one lincRNA and two circRNAs that could sustain high expression of the core transcription factor Oct4 by sponging miRNAs during reprogramming, and they ultimately improved the reprogramming efficiency and quality of iPSCs. The knockdown of each of the lncRNA could significantly down-regulate Oct4 expression, but the formation of AP-positive colonies was not affected, indicating the lncRNAs as miRNA sponges can mutually compensate and iPSCs reprogramming can tolerate a certain downregulation of exogenous Oct4 expression. Furthermore, considering that the formation of AP-positive colonies was not significantly changed by knockdown or overexpression of one single lncRNA, we suggest that the combined effect of the lncRNAs as miRNA sponges, rather than a single lncRNA, plays the central role on iPSCs reprogramming. During reprogramming, the activation of endogenous pluripotent genes happens at late stages (27,36). So, here, the ceRNA mechanism is involved in the regulation of exogenous Oct4 expression. Whether this mechanism regulates endogenous pluripotent gene expression during reprogramming During reprogramming, the miRNAs targeting Oct4 mRNA from MEFs are highly expressed, however, the expression of exogenous Oct4 is not suppressed. Linc/circRNAs, which have complementary binding sites with the miRNAs targeting on Oct4 are also highly expressed and can counteract the activities of the miRNAs as sponges. The knockdown of the linc/circRNAs results in down-regulation of Oct4 expression, the imprinting defect at the Dlk1-Dio3 locus, decrease of reprogramming efficiency, and low-grade chimera forming iPSCs, and those can be rescued by Oct4-addback. In summary, the linc/circRNAs can improve reprogramming efficiency and quality of iPSCs by sponging miRNAs targeting on Oct4 during reprogramming.

LincRNAs and circRNAs sponge miRNAs during reprogramming
remains to be further determined. Together with previous studies (11), we suggest that the ceRNA network plays an important role on establishment and maintenance of pluripotency during reprogramming and in pluripotent stem cells.
In summary, our results point out a high prevalence of miRNA-mediated interactions between their mRNA targets and lncRNAs, proposing that this mechanism of ceRNAs indeed functions on regulation of gene expression, particularly in the context of cell-fate transitions.

Animals
Mice were housed and prepared according to the protocol approved by the Institutional Animal Care and Use Committee of Northeast Agriculture University (protocol number IACUC-02-005) and the IACUC of Weill Cornell Medical College (protocol number 2014-0061).

High throughput RNA-seq
Total RNA was isolated from the cells at days 0, 3, and 6 post-OSKM induction and the fully reprogrammed iPSCs using TRIzol RNA extraction reagent (Thermo Fisher Scientific), following the manufacturer's protocol. The concentration and purity of total RNA were assessed on Nanodrop 2000 (Thermo Fisher Scientific), and the integrity of total RNA was evaluated with 2% agarose gel electrophoresis and a RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA). The cDNA libraries of long chain RNAs (mRNA/lncRNA/circRNA) and small RNAs were generated from 10 and 2 g of total RNA, respectively. rRNA was removed using a Epicenter Ribo-zero rRNA Kit (Epicenter). Sequencing libraries were generated using the NEBNext Ultra TM

RNA-seq data processing and analysis
The qualified reads were aligned with Tophat 2 into the mouse reference genome by default parameters. The reads on the extracted alignment were constructed with Cufflinks to construct transcripts. The transcripts of all samples were combined and reconstructed into a large transcript file (merged.gtf) using Cuffmerge. Then the expression of genes in each sample were quantified to FPKM. For differential expression, we first counted the overlap of reads with genes by htseq-count. Next, we compared the two groups using default parameters in R package DESeq. A gene was considered significant if the Benjamini and Hochberg-adjusted p value (Padj) was less than 5% and the fold-change was greater than 2.

Quantification of miRNA abundance and prediction of miRNA response elements
Clean sequencing reads were obtained by removing from the raw data reads containing more than one low quality (Q-value Յ20) base, with 5Ј primer contaminants, without the 3Ј primer, without the insert, with poly(A), and shorter than 18 nucleotides. Clean reads were aligned with the reference genome (ftp:// ftp.ensembl.org/pub/release-92/fasta/mus_musculus/dna/ Mus_musculus.GRCm38.dna.toplevel.fa.gz), the Rfam database 13.0 (37), and the RepBase (38) to identify and remove rRNA, scRNA, snoRNA, snRNA, tRNA, repeat sequences, and fragments from mRNA degradation. The remaining reads were searched against miRBase 22.0 (39) to identify known miRNAs in mouse. The possible MREs of known miRNAs were predicted by RNA22 (32), an miRNA target prediction tool, with default parameter setting and the binding energy cutoff was less than or equal to Ϫ20 Kcal/mol. The numbers of MREs per kb length of each upregulated or down-regulated gene and the average numbers of MREs in each group were calculated. The miRNA expression level was calculated and normalized using RPM (reads per million). We identified significant differentially expressed miRNAs with a false discovery rate Ͻ0.05 and fold-change Ն2 threshold by DEGseq (40).

Identification of linc/circRNAs and expression analysis
Pre-processing of raw sequencing data were performed using the FASTX-Toolkit with default parameters by removing low quality reads (more than 20% of the bases qualities are lower than 10), reads with adaptors, and reads with unknown bases (n bases more than 5%). Clean reads were mapped to the mouse reference genome (ftp://ftp.ensembl.org/pub/release-92/fasta/ mus_musculus/dna/Mus_musculus.GRCm38.dna.toplevel.fa.gz). Unmapped reads were collected and the software CIRI (41) and find_circ (17) were used to identify circRNAs with default options, and circBase IDs were used to indicate known circRNAs. Counts of identified circRNA reads were normalized by the RPB (junction reads per billion mapped reads) method. For lincRNAs, the reconstruction and identification of transcripts were performed by mapped reads with StringTie (42) and cuffcompare (43). Known lincRNAs were acquired by gene annotations for comparing the assembled transcripts with NCBI, Ensembl, and UCSC mouse known genes. The expression level of lincRNA was calculated using FPKM (fragments per kilobase of transcript per million mapped reads) with the software RSEM (44). Transcripts with a false discovery rate Ͻ0.05 and fold-change Ն2 were then identified as significant differentially expressed lincRNA or circRNA using DEGSeq (40).

Experimental validation of miRNA activity and interaction between mRNAs and linc/circRNAs
The miRNAs mimics and inhibitors, and linc/circRNAs siRNAs were synthesized by Exiqon, and transfected in triplicates into MEFs using Lipofectamine RNAi MAX Reagent (Invitrogen) according to the manufacturer's guidelines. rMEFs were collected at day 3 post-OSKM induction, RNAs were extracted and checked by qRT-PCR for gene expressions. RIP was performed using a Magnetic RNA-Protein Pull-down Kit (Pierce, number 20164) according to the manufacturer's protocol. Briefly, the 3Ј end-biotinylated miRNAs mimics were transfected into MEFs at a final concentration of 20 nM for 1 day. The biotin-coupled RNA complex was pull downed by incubating the cell lysates with streptavidin-coated magnetic beads. The abundance of linc/circRNAs in bound fractions was evaluated by qRT-PCR analysis normalized to Oct4 mRNA. Gapdh was analyzed as negative control. The TetO-FUW-Oct4 plasmid (Addgene plasmid 20323) was used for Oct4-addback, and 10,000 cells infected with a dose of 5 ϫ 10 8 cfu of virus for 3 h.

Methylation analysis of the Dlk1-Dio3 locus
Genomic DNA isolation from cells and bisulfite pyrosequencing were performed as follows. Briefly, genomic DNA of iPSCs was extracted with the Gentra Pure gene kit (Qiagen) according to the manufacturer's instructions and quantified using a Nano Drop 2000 spectrophotometer. Bisulfite conversion of 500 ng of genomic DNA per sample was performed with the EpiTect Bisulfite Kit (Qiagen) according to the manufacturer's specifications. Quantification of DNA methylation was carried out by PCR of bisulfite-converted DNA and pyrosequencing. PCR and sequencing primers for bisulfite pyrosequencing were designed using Pyrosequencing Assay Design Software 2.0 (Qiagen). For pyrosequencing, a PyroMark Q96 ID instrument (Qiagen) and PyroMark Gold Q96 reagents (Qiagen) were used. Data were analyzed using the PyroMark CpG Software 1.0.11 (Qiagen).

Statistical analysis
All data are presented as mean Ϯ S.D. Differences between groups were tested for statistical significance using Student's t test or 2 test. Statistical significance was set at p Ͻ 0.01.