Common microRNA–mRNA interactions exist among distinct porcine iPSC lines independent of their metastable pluripotent states

Previous evidences have proved that porcine-induced pluripotent stem cells (piPSCs) could be induced to distinctive metastable pluripotent states. This raises the issue of whether there is a common transcriptomic profile existing among the piPSC lines at distinctive state. In this study, we performed conjoint analysis of small RNA-seq and mRNA-seq for three piPSC lines which represent LIF dependence, FGF2 dependence and LFB2i dependence, respectively. Interestingly, we found there are 16 common microRNAs which potentially target 13 common mRNAs among the three piPSC lines. Dual-luciferase reporter assay validated that miR-370, one of the 16 common microRNAs, could directly target the 3′UTR of LIN28A. When the differentiation occurred, miR-370 could be activated in piPSCs and switched off the expression of LIN28A. Ectopic expression of miR-370 in piPSCs could reduce LIN28A expression, decrease the alkaline phosphatase activity, slow down the proliferation, and further cause the downregulation of downstream pluripotent genes (OCT4, SOX2, NANOG, SALL4 and ESRRB) and upregulation of differentiation relevant genes (SOX9, JARID2 and JMJD4). Moreover, these phenotypes caused by miR-370 could be rescued by overexpressing LIN28A. Collectively, our findings suggest that a set of common miRNA–mRNA interactions exist among the distinct piPSC lines, which orchestrate the self-renewal and differentiation of piPSCs independent of their metastable pluripotent states.

MicroRNAs (miRNAs) are a group of small noncoding RNAs that consist of about 20-21 nucleotides and function in post-transcriptional regulation of gene expression. Several miRNA clusters, including miR-290-295/miR-371-373, miR-302 and miR-17-92b, are defined as the embryonic stem cell cycle-specific (ESCC) miRNAs. 1 Recent evidences have showed that these ESCC miRNAs have an important role in maintaining pluripotency and cellular reprogramming. 2 The miR-302/367 cluster was reported to regulate both cell cycle and apoptosis in human embryonic stem cells (hESCs). 3 In addition, overexpression of miR-302a, miR-302b and miR-200c could improve the efficiency of porcine-induced pluripotent stem cells (piPSCs) induction. 4 To further study the post-transcriptional regulation mechanism of global miRNAs on pluripotency, small RNA sequencing has been used extensively to reveal transcription profiles. The differentially expressed miRNAs have been documented between human-induced pluripotent stem cells (hiPSCs) and hESCs, 5 in which miR-371-3 are highly expressed in hESCs versus hiPSCs. In contrast, miR-181a, miR-199b-3p and miR-214 are hiPSC-specific miRNAs, which are lowly expressed in hESCs. In addition, comparison of global miRNA profiles also showed differences between naive state and primed state on mice. 6 Although many attempts have been tried to get piPSCs similar to mouse naïve pluripotent state, it is still challenging to achieve authentic piPSCs that can form germline-transmission-chimera with high efficiency. However, distinctive piPSC lines have been made, which featured by LIF dependence, 7 FGF2 dependence 8 and LFB2i dependence. 9 The LIF-dependent piPSCs were generated by transducing four Yamanaka factors into porcine embryonic fibroblasts (PEFs) and propagated in the presence of leukemia inhibitor factor (LIF) and forskolin. 7 By contrast, the FGF2-dependent piPSCs were induced by ectopically expressing four Yamanaka factors plus LIN28 and c-MYC into porcine mesenchymal stem cells and maintained in the presence of fibroblast growth factor 2 (FGF2). 8 Our previous study showed that LFB2idependent piPSCs could be made by overexpressing four Yamanaka factors into PEFs and maintained with addition of three growth factors (LIF, FGF2 and BMP4) and two inhibitors (2i: CHIR99021 and SB431542). 9 Although these three types of piPSC lines were generated under diverse cellular context and by different reprogramming strategies, they represent different metastable pluripotent states. Recently, the miRNA profiles of piPSCs derived from culture conditions of both hESCs and mouse embryonic stem cells (mESCs) were investigated and showed the significant differences between the two culture conditions. 10 However, it is unclear whether there is a common miRNA profile existing among the distinct piPSC lines.
In this study, we performed conjoint analysis of small RNAseq and mRNA-seq for three metastable piPSC lines with dependence of distinctive signaling pathways. Interestingly, we found there were 16 common miRNAs that potentially target 13 common mRNAs among the three piPSC lines. MiR-370, one of the 16 miRNAs, was demonstrated to target the 3′UTR of LIN28A. Ectopic expression of miR-370 could reduce the expression level of LIN28A in piPSCs. Even more, miR-370 could be activated and switched off the expression of LIN28A upon piPSCs differentiation. We also found that ectopic expression of miR-370 in piPSCs could decrease the alkaline phosphatase activity, slow down the cellular proliferation, and further cause the downregulation of downstream pluripotent genes (OCT4, SOX2, NANOG, SALL4 and ESRRB) and upregulation of differentiation relevant genes (SOX9, JARID2 and JMJD4). Of note, these defects caused by miR-370 could be rescued by overexpressing LIN28A.Our study uncovered that a set of common miRNA-mRNA interactions exist among distinct piPSC lines independent of their metastable pluripotent states. The orchestration of miRNA-mRNA interactions is supposed to regulate the self-renewal and differentiation of piPSCs.

Results
Global miRNA expression pattern of distinct piPSC lines. To determine whether there is a common miRNA profile existing in different types of piPSCs, three piPSC lines that represent LIF dependence, FGF2 dependence and LFB2i dependence were selected for small RNA-seq analysis. The reprogramming context of the three piPSC lines and sequencing information is summarized in Figure 1a. The small RNA-seq of PEFs was used as a somatic cell control. Small RNA-seq showed the number of miRNAs and its proportion in total small RNAs for LIF-dependent piPSCs (piPS-L), FGF2-dependent piPSCs (piPS-F), LFB2idependent piPSCs (piPS-LF) and PEFs is 11144 (80.85%), 6838 (69.92%), 10146 (44.92%) and 7661 (75.03%), respectively ( Figure 1b). To annotate the acquired pig miRNAs, Blast analysis was performed based on pig miRBase as well as other animal species' miRbases, including human, mouse, cattle and rat. 11 Results showed that only 21.2 to 34.6% of porcine miRNAs was annotated. In contrast, most sequenced miRNAs (~78%) have no annotation due to the limited information of pig miRBase (Figure 1c). Besides, we compared the differentially expressed miRNAs between each line of piPSCs and PEFs. The number of differentially expressed miRNAs (FDRo0.05) in piPS-L, piPS-F and piPS-LF is 897, 1068 and 306, respectively (Figure 1d; Supplementary Table 1). Then, we carried out Venn analysis of the differentially expressed miRNAs and revealed that 470 miRNAs were specific for piPS-L, 651 for piPS-F, and 67 for piPS-LF. We also found that 220 miRNAs were shared by piPS-L and piPS-F, 32 for piPS-F and piPS-LF, and 42 for piPS-L and piPS-LF. Importantly, there were 165 miRNAs which were commonly expressed among the three piPSC lines (Figure 1e; Supplementary Table 2). Collectively, the result of small RNA-seq revealed that a common miRNA profile existed among distinct piPSC lines, although piPSC lines differed in their signaling dependence and original context.
Analysis of common miRNAs revealed by small RNAseq. The heatmap of 165 common miRNAs showed that piPSC lines were clustered relatively closer versus PEFs. Three miRNAs (miR-371, chr8_17989_mature and chr15_31797_mature) were highly expressed in piPSC lines, but lowly expressed in PEFs. By contrast, the other 162 miRNAs showed the lower expression level in piPSC lines versus PEFs (Figure 2a). The three high-expressed miRNAs in piPSCs were predicted to target four mRNAs, while 162 lowly expressed miRNAs in piPSCs were predicted to target 682 mRNAs (Figure 2b; Supplementary Table 3). We then performed GO and KEGG analysis of the predicted 682 mRNAs. The enriched biological processes of GO analysis (Po0.05) included in utero embryonic development, positive regulation of peptidyl-tyrosine phosphorylation, actin cytoskeleton reorganization, axonal fasciculation, autophagy, positive regulation of Notch signaling pathway, glutathione metabolic process and retrograde transport, and endosome to Golgi (Figure 2c; Supplementary Table 4). The significantly enriched KEGG pathways (Po0.05) included metabolic pathways, histidine metabolism, MAPK signaling pathway, arachidonic acid metabolism, taurine and hypotaurine metabolism, glycerophospholipid metabolism and fatty acid metabolism (Figure 2d; Supplementary Table 4). These results indicated that a set of miRNAs were commonly expressed with low level in piPSC lines and they were supposed to target on genes involved in many important biological processes and pathways.  Table 6). We then focused on the common mRNAs that were divisible into two categories: 391 highly expressed mRNAs and 979 lowly expressed mRNAs (Figure 3b; Supplementary Table 7). GO analysis and KEGG analysis were performed for these common mRNAs (Figures 3c and d; Supplementary Table 8). The common high-expression mRNAs are enriched in typical GO terms and pathways of stem cells, such as 'stem cell differentiation' and 'tight junction'. In contrast, the common low-expression mRNAs are mostly enriched in GO terms and pathways relevant to extracellular matrix. This observation uncovered fundamental transcriptome differences between piPSCs and somatic cells.
Conjoint analysis of small RNA-seq and mRNA-seq for distinct piPSC lines. To increase the accuracy of predicted targets from the common miRNAs for distinct piPSC lines, conjoint analysis of small RNA-seq and mRNA-seq was performed and validated by wet-lab experiments (Figure 4a).  Table 9). These common mRNAs include CHGA, EPB414A, OTX2, TARSL2, LOC100521376, LIN28A, CYP2D25, CAMSAP3, SASH3, RGS4, RNF207 and RAB33A. Interestingly, two well-known pluripotent genes OTX2 and LIN28A were found among the list. The interaction network of the common miRNAs and these common 13 mRNAs was constructed ( Figure 4c). Interestingly, miR-206 was found to target OTX2 while miR-370 and chr15_31863_mature (miR-31863) were found to target LIN28A. The negative correlation of miR-370 and miR-31863 with LIN28A was confirmed in another cell line of piPSCs (Dox-piPSCs) described in this lab 12 and PEFs by qRT-PCR, showing that expression level of miR-370 and miR-31863 were high and LIN28A was low in PEFs, conversely the expression level of miR-370 and miR-31863 were low and LIN28A was high in piPSCs ( Figure 4d). The MiR-370 targets the 3′UTR of LIN28A mRNA. The following study focused on the interaction of miR-370 and miR-31863 with LIN28A in piPSCs. Because miR-31863 is a novel miRNA found in pig, its secondary structure was predicted by using RNAFold and we found pre-miR-31863 could form a stem loop structure similar to pre-miR-370 ( Figure 5a). The analysis showed that LIN28A 3′UTR retains two potential sites for miRNA targeting. One is for miR-31863 and the other one is for miR-370 ( Figure 5b). To determine the actual interaction, a reporter vector that contains the showed that mimics of miR-370 (mimics-370) significantly downregulated luciferase activity, but mimics-31863 had no effect (Figure 5d). Then we focused on mimics-370 and examined its dose-dependent effect on LIN28A 3′UTR reporter. As the concentration of mimics-370 increased, the luciferase activity decreased gradually with dose-dependent course (Figure 5e). We also investigated the downregulation effect of luciferase activity by mimics-370 over the timecourse. However, no significant difference was observed during the treatment time from 24 h to 48 h (Figure 5f). To determine the specificity of mimics-370 targeting on LIN28A, we constructed a reporter vector that contained the mutated sequence of 3′UTR of LIN28A (MUT-UTR) (Figure 5g). The cotransfection result showed that the MUT-UTR completely abolished the mimics-370 effect when compared with WT-UTR ( Figure 5g). These results suggest that miR-370 can directly target the 3′UTR of LIN28A mRNA.
MiR-370 represses LIN28A expression. To validate the regulation of LIN28A expression by miR-370, we investigated expression changes of miR-370 and LIN28A in Dox-piPSCs, a tetracycline operator (TetO)-inducible piPSCs line reported in this laboratory previously. 12 The Dox-piPSCs were maintained in the pluripotent state with an addition of doxycycline, but started to differentiate soon after withdrawal of doxycycline ( Figure 6a). In the differentiated Dox-piPSCs, miR-370 was significantly increased and the expression of LIN28A was significantly decreased (Figure 6b). This finding intrigued us to investigate which transcription factors may be potentially involved in activating miR-370 when the Dox-piPSCs initiated differentiation. We analyzed the promoter of Pri-miR-370 in the DLK1-MEG3 locus 13 and found that the genomic DNA sequences of Pri-miR-370 promoter region were highly conserved among pig, human and mouse ( Figure 6c). We further predicted transcription factor binding sites in the Pri-miR-370 promoter of pig, human and mouse, respectively, by searching JASPAR database.  Table 12). Together these 18 transcription factors are suggested to potentially activate miR-370 when the piPSCs initiate differentiation.
To validate the repression of LIN28A expression by miR--370, we next transduced the piPS-LF with miR-370 lentiviruses. The piPS-LF colonies with miR-370 overexpression became morphologically incompact even cultured under LFB2i condition 9 (Figure 6d). The endogenous LIN28A was significantly reduced on both the mRNA level and protein level in piPS-LF upon overexpressing miR-370 (Figures 6e and f). To test whether the interaction between LIN28A and miR-370 exists universally in other piPSCs, we further investigated the effects of miR-370 on Dox-piPSCs. The Dox-piPSCs transduced with miR-370 lentiviruses showed similar changes which were observed in piPS-LF (Figures 6g-i). Then we extended our investigation of interaction between Lin28a and miR-370 in mESCs line J1. The 3′UTR of mouse Lin28a mRNA also contains a targeting site for miR-370 (Supplementary Figure 2). So we asked whether overexpression of miR-370 in mESCs could also lead to similar phenotype as observed in piPSCs. Unexpectedly, the mESCs with miR-370 overexpression showed no obvious morphology change (Figure 6j), although the expression of Lin28a was significantly decreased (Figures 6k and l). Together the results above indicate that miR-370 can conservatively repress LIN28A expression in diverse pluripotent stem cells instead of a species-specific manner.
MiR-370 affects self-renewal and differentiation of piPSCs by LIN28A dependence. To further explore the effect of miR-370 on self-renewal and differentiation of piPSCs, we firstly investigated the impact of miR-370 on alkaline phosphatase (AP) activity. Interestingly, we found that miR-370 overexpression decreased the AP activity in Dox-piPSCs. However, LIN28A overexpression rescued the weakened AP activity in the Dox-piPSCs caused by miR-370 (Figure 7a).
Then we examined the impact of miR-370 on cellular proliferation. We noted that miR-370 overexpression slowed down the proliferation rate of Dox-piPSCs. Similarly, this proliferation defect caused by miR-370 could be rescued by LIN28A overexpression in Dox-piPSCs ( Figure 7b).
Finally, we observed the impact of miR-370 on the expression of pluripotency and differentiation relevant genes. We found that miR-370 overexpression reduced the expression level of pluripotent genes (LIN28A, OCT4, SOX2, NANOG, SALL4 and ESRRB) in Dox-piPSCs. By contrast, miR-370 overexpression in Dox-piPSCs upregulated the expression level of differentiation relevant genes (SOX9, JARID2 and JMJD4) (Figure 7c). Of note, LIN28A overexpression could rescue the downregulation of pluripotent genes and upregulation of differentiation genes caused by miR-370 in Dox-piPSCs (Figure 7c). We also investigated the impact of miR-370 on the expression of pluripotent genes in piPS-LF and mESCs. Similar phenotype was observed   Figure 3).
Taken together, these results indicate that miR-370 can regulate the self-renewal and differentiation of piPSCs by LIN28A dependence.

Discussion
In this study, we proved that a set of common microRNA-mRNA interactions exists among distinct piPSC lines independent of their metastable pluripotent states. Furthermore, we validated the interaction of miR-370 with LIN28A, and demonstrated that miR-370 can target the 3′UTR of LIN28A and reduce the LIN28A expression in piPSCs. We propose the working model of interaction between miR-370 and LIN28A in piPSCs (Figure 7d). In the pluripotent state, LIN28A is highly expressed while miR-370 is lowly expressed. LIN28A protein is speculated to inhibit the maturation of miR-370 precursor. 14 When the differentiation initiates, miR-370 is highly activated and inhibits LIN28A expression  MiRNAs represent an important mechanism of posttranscriptional regulation and are highly conserved between invertebrates and vertebrates. 15 We revealed there are 165 common miRNAs that are differentially expressed between distinct piPSC lines and somatic cells. Interestingly, only three out of the 165 common miRNAs are highly expressed in piPSCs and others are highly expressed in somatic cells. The three miRNAs include miR-371 and other two predicted miRNAs (chr8_17989_mature and chr15_31797_mature). This unique miRNA signature of piPSCs presented an obvious contrast against the specific miRNA signature of mouse and human PSCs. Substantial evidence has shown that miRNA clusters miR-302, miR-290 and miR-371 are expressed exclusively in both mouse and human pluripotent stem cells (PSCs). 16 Taken the findings together, miR-371 seems to be more conserved than miR-302 and miR-290 in PSC across diverse species. Recent evidence suggested that miR-371 could promote pluripotency by regulating glycolytic metabolism via the Mbd2-Myc signaling pathway in mouse and human ESCs. 17 Thus, it will be interesting to investigate if the analogous role of miR-371 exists in piPSCs in the future.
We also uncovered that there are 1416 common mRNAs differentially expressed between distinct piPSC lines and somatic cells. Out of the 1416 common mRNAs, many genes are commonly expressed in somatic cells with high levels and enriched in biological processes of extracellular matrix organization and cell adhesion. This feature of gene expression pattern reflects the fundamental differences on extracellular matrix between PSC and somatic cells. Extracellular matrix organization and cell adhesion of somatic cells are regarded as barriers to achieve fully reprogramming. 18 Moreover, conjoint analysis of small RNA-seq and mRNAseq was applied and revealed that sixteen miRNAs were commonly expressed lower in piPSC lines than that in somatic cells, and these miRNAs potentially targeted thirteen common mRNAs. These miRNA-mRNA interactions are common in distinct piPSC lines independent of their metastable pluripotent states. Interestingly, interactions of miR-206 on OTX2 and miR-370 on LIN28A were found among the common interactions. Otx2, once as a transcription factor involved in brain development, 19 was found recently to antagonize naïve pluripotency and promote the naïve-to-primed transition. 20 We previously demonstrated that overexpression of OTX2 decreased the percentage of alkaline phosphatase-positive piPSC colonies and downregulated the expressions of NANOG and OCT4. 21 The common interaction of miR-206 on OTX2 found in this study suggests the potential use of miR-206 to induce piPSCs with naïve pluripotency by downregulating OTX2.
LIN28A, as an RNA-binding protein, has important roles in tumorigenesis, embryonic development and pluripotency regulation. 22 Overexpression of LIN28A and LIN28B could increase iPSCs derivation efficiency, while double knockout of LIN28A/B reduced the reprogramming efficiency and caused mouse iPSCs trapped in a more naïve state. 23 In addition, mESCs at naïve state showed the low expression of Lin28a, whereas mouse epiblast stem cells at primed state express Lin28a with high level. 24 Lin28a was further demonstrated to facilitate the naive-to-primed transition in mouse PSCs. 23 All these evidences indicate that LIN28A can regulate PSCs exiting from naïve pluripotency as a key factor. Our findings showed that LIN28A was conservatively expressed in distinct piPSC lines, which indicated the primed pluripotent feature of these cell lines and might explain the limited chimera potential for these piPSCs. [7][8][9] Moreover, we demonstrated that miR-370 could target on the 3′UTR of LIN28A and reduce the expression of LIN28A in piPSCs. Therefore, ectopic expression of miR-370 may logically push piPSCs to acquire naive pluripotency. Unexpectedly, we found ectopic expression of miR-370 could not make piPSCs acquire naïve pluripotency. In contrast, miR-370 overexpression led to the downregulation of classical pluripotent genes (including naïve pluripotency gene NANOG) in piPSCs. Taken these findings together, LIN28A seems to have different roles in PSC across species. Additionally, we also demonstrated that miR-370 was activated during the process of piPSCs differentiation and decreased the expression level of LIN28A, which suggests the role of miR-370 as an orchestrator for regulating self-renewal and differentiation of piPSCs. The similar interaction of miR-370 on LIN28A was once reported in human hepatocellular carcinoma and demonstrated that LIN28A protein could blocked the biogenesis of miR-370 by an inverse feedback. 14 It will be interesting to investigate this inverse feedback effect of LIN28A on miR-370 in piPSCs (Figure 7d) in the future study.
Collectively, a set of common microRNA-mRNA interactions was uncovered among distinct piPSC lines independent of their metastable pluripotent states, which will provide insights into the mechanism of reprogramming and pluripotency regulation on pig.
Small RNA sequencing. For small RNA sequencing, the RNA lysates of piPS-L and piPS-F was provided by Dr. Hanazono's lab and Dr. West's lab, respectively. The RNA lysates of piPS-LF was prepared by our lab. Three microgram RNA per sample was used as input material for the small RNA library. Sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, MA, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, NEB 3′ SR Adaptor was ligated to 3′ ends of miRNA, small interfering RNA (siRNA) and piwi-interacting RNA (piRNA) directly and specifically. After the 3′ ligation reaction, the SR RT Primer hybridized to the excess of 3′ SR Adaptor (that remained free after the 3′ ligation reaction) and transformed the single-stranded DNA adaptor into a double-stranded DNA molecule. The 5'end adapter was ligated to 5′ ends of miRNAs, siRNA and piRNA. Then first-strand cDNA was synthesized using M-MuLV Reverse Transcriptase (RNase H -). PCR amplification was performed using LongAmp Taq 2X Master Mix, SR Primer for illumina and index (X) primer. PCR products were purified on 8% polyacrylamide gel (100 V, 80 min). DNA fragments corresponding to 140-160 bp (the length of small noncoding RNA plus the 3′ and 5′ adaptors) were recovered and dissolved in 8 μl elution buffer. The library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia, San Diego, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500/2000 platform and 50 bp single-end reads were generated. Small RNA-seq data of piPS-L, piPS-F, piPS-LF and PEFs were deposited in EMBL-EBI database (http://www.ebi.ac.uk/) under accession number E-MTAB-2631. mRNA sequencing. For mRNA sequencing, a total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First-Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/ polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with the hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (New England Biolabs) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 95°C 5 min before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified using AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on Illumina Hiseq 2000 platform and 100 bp pairedend reads were generated. The sequenced data have been deposited with the European Bioinformatics Institute (www.ebi.ac.uk/). The mRNA-seq data of piPS-L, piPS-LF and PEFs were deposited in EMBL-EBI database (http://www.ebi.ac.uk/) under accession number E-MTAB-2634. mRNA-seq data of piPS-F were downloaded from GEO under the accession number GSM881389.
Data analysis. After removal of adaptors, low quality tags, and contaminants from the sequenced tags, clean reads were annotated. Index of the reference genome was built using Bowtie v2.0.6. 25 Paired-end clean reads were aligned to the reference genome using TopHat v2.0.9. 26 Normalization was performed before further analysis. Hierarchical clustering was done by using gplots package. The miRNA profiles were clustered using average linkage clustering with Euclidean distances, treating samples independently each other. Differentially expressed miRNAs between any two samples were calculated by EBSeq. 27 Target genes of miRNAs were predicted using RNAhybrid. 28 Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed genes were performed using the Database for Annotation and Visualization and Integrated Discovery (DAVID); 29,30 a hypergeometric test with the Benjamini and Hochberg false discovery rate (FDR) was performed using the default parameters to adjust the P-value.
Secondary structure of miRNA was predicted by RNAFold. 31 The promoter region of miR-370 was predicted by the Web Promoter Scan Service (http://www-bimas.cit. nih.gov/molbio/proscan/). The potential transcription factors binding to the miR-370 promoter was predicted by the JASPAR. 32 Gene cloning and vector constructs. Genomic DNA (gDNA) was isolated from PEFs and its concentration and purity was examined using Nanodrop (Thermo Scientific, San Jose, CA, USA). The 3′UTR of pig LIN28A were amplified by PCR using pig gDNA and the mutated UTR was amplified by Overlap PCR. The cloned fragments were cut with XhoI and NotI, and ligated into pSICHECK2 vector (Promega, Madison, USA). The synthesized pre-mir-370 (80 bp) was cloned into pLL3.7 basic (Addgene, Cambridge, MA, USA) with HpaI and XhoI. The CDS of porcine LIN28A and mouse Lin28a were amplified by PCR using mRNA isolated from DOX-piPSCs and J1. The cloned fragments were cut by EcoRI and BamHI, and ligated into pCD513B-1 vector (System Biosciences, Mountain View, CA, USA). The synthesized mouse pri-miR-370 (537 bp) fragment was also cloned into pCD513B-1 vector. All the primers used for vector construction are listed in Supplementary Table 13.
Cell transfection and luciferase assay. All miRNA mimics were ordered from GenePharma. Cells were cotransfected by miRNA mimics reporter plasmids using X-tremeGENE siRNA Transfection Reagent (Roche, Mannheim, Germany). Cells were harvested 36 h after transfection and lysed for 15 min at room temperature. The firefly and renilla luciferase activities were analyzed using Dual-luciferase reporter assay system (Promega, Wisconsin, USA) on a Hamamatsu BHP9504 Luminometer (Hamamatsu, Peking, China). Firefly luciferase activity was normalized to the renilla luciferase activity.
Lentivirus package and infection. HEK 293T cells were seeded in a 6-well plate the day before transfection. After 24 h culture, cells were transfected with lentivirus vector (pLL3.7-miR-370, pCDH-miR-370, pCDH-LIN28A (pig), or pCDH-Lin28a (mouse)) and its package plasmids (pSPAX2, pVSVG) using Lipofectamine 2000 (Invitrogen), according to the manufacturer's instruction. At 48 h after transfection, the supernatant with lentivirus particles was collected and used for infecting Dox-piPSCs, piPS-LF or J1 in the presence of 4 ug/ml polybrene (Sigma-Aldrich). Two rounds of infection were conducted. The infected piPSCs colonies with GFP or puromycin resistance were picked up and expanded for further analysis.
Quantitative RT-PCR. Total RNAs from Dox-piPSCs were reverse transcribed to cDNAs using Revert Aid First-Strand cDNA Synthesis Kit (Thermo Scientific) and Mir-X™ miRNA First-Strand Synthesis and SYBR qRT-PCR (Clontech, Palo Alto, CA, USA) according to the manufacturer's instruction. For the quantitative RT-PCR, reactions were performed in triplicate using SYBR Green PCR Master Mix (TransGen Biotech, Peking, China), and detected with the CFX96 real-time PCR system (Bio-Rad, Hercules, USA). Relative expression levels of target genes and deferentially expressed miRNAs were normalized to either U6 (for miRNAs) or β-Actin expression (for mRNAs). The relative expression levels were calculated using 2 − ΔΔCt .
Alkaline phosphatase staining. The alkaline phosphatase (AP) activity of piPSCs was determined by AST Fast Red TR (Sigma-Aldrich) and a-Naphthol AS-MX Phosphate (Sigma-Aldrich) according to the manufacturer's instruction. Briefly, cells were washed twice with ice-cold PBS, fixed with 4% paraformaldehyde (pH 7.4) for 15 min at room temperature, and followed by washing three times with icecold PBS. Then the cells were then incubated at room temperature with dye solution containing AST Fast Red TR (1.0 mg/ml), a-Naphthol AS-MX (0.4 mg/ml) in 0.1 M Tris Buffer. After 10-20 min incubation, the cells were washed three times with PBS and the images were documented with a Nikon microscope.
Cell growth curve. Growth curve of Dox-piPSCs was determined by cell counting for five days. The cells were plated with 4 × 10 4 cells/well in 12-well plate and harvested for cell counting every 24 h. The cells in triplicate wells were counted and averaged for each treatment at every time point.
Statistics analysis. Calculated data are presented as mean ± S.D. Student t-test was used to determine the differences between two groups and one-factor ANOVA with Bonferroni post hoc test was used to determine the differences among three groups in this study. All statistical analyses were done with GraphPad Prism 5.0 (GraphPad, La Jolla, CA, USA).