Oct4-mediated reprogramming induces embryonic-like microRNA expression signatures in human fibroblasts

Oct4-mediated reprogramming has recently become a novel tool for the generation of various cell types from differentiated somatic cells. Although molecular mechanisms underlying this process are unknown, it is well documented that cells over-expressing Oct4 undergo transition from differentiated state into plastic state. This transition is associated with the acquisition of stem cells properties leading to epigenetically “open” state that is permissive to cell fate switch upon external stimuli. In order to contribute to our understanding of molecular mechanisms driving this process, we characterised human fibroblasts over-expressing Oct4 and performed comprehensive small-RNAseq analysis. Our analyses revealed new interesting aspects of Oct4-mediated cell plasticity induction. Cells over-expressing Oct4 lose their cell identity demonstrated by down-regulation of fibroblast-specific genes and up-regulation of epithelial genes. Interestingly, this process is associated with microRNA expression profile that is similar to microRNA profiles typically found in pluripotent stem cells. We also provide extensive network of microRNA families and clusters allowing us to precisely determine the miRNAome associated with the acquisition of Oct4-induced transient plastic state. Our data expands current knowledge of microRNA and their implications in cell fate alterations and contributing to understanding molecular mechanisms underlying it.


Results
Oct4 over-expression down-regulates mesenchymal and fibroblast-specific genes and up-regulates epithelial genes. We transduced human dermal fibroblasts (hDFs) with lentiviral vectors for Oct4 over-expression and GFP for control experiments. We chose hDFs, because they are readily obtainable from patient skin biopsies, easily cultured in vitro and typically used as a starting cell line in reprogramming experiments. To determine the nature of Oct4-induced plasticity, hDFs were cultured in hDFs media in the absence of lineage-inducing growth factors that have been used in previous studies [4][5][6] , as the presence of these factors would bring bias to our analysis. Cells were harvested 6 days upon transduction and Puromycin selection. We chose this time-point, because 6 days provides enough time for antibiotic selection to yield homogenous population of cells expressing Oct4.
Transduced hDFs over-expressed Oct4, showed a dramatic change of morphology shortly upon Oct4 over-expression with transition of long-spindled fibroblast morphology to short-spindled cell shape (Fig. 1a,b), and maintained this altered morphology for at least 30 days (Fig. S1). In order to further analyse molecular mechanisms underlying changed morphology, we aimed to assess the expression of mesenchymal and epithelial genes. Western blot analysis revealed down-regulation of mesenchymal and fibroblast markers such as Slug, N-cadherin, Vimentin and up-regulation of epithelial marker ZO-1 (Figs 1c, S2a). Interestingly, we also detected up-regulation of Snail (SNAI1), which is classically viewed as a mesenchymal gene, however, we and others have previously shown that Snail is up-regulated during early stages of reprogramming and enhances efficiency of iPSCs generation 25,30 . RT-qPCR analysis showed down-regulated mesenchymal genes COL1A1, SNAI2 (SNAI2 not significantly) and not significantly up-regulated epithelial genes CDH1, EPCAM, and CRB3 upon Oct4 over-expression (Fig. 1d).
Given the observed change of cell morphology and altered mesenchymal/epithelial gene expression, we sought to further investigate, if Oct4 over-expression affects cell migration. Scratch-wound assay showed that control (GFP+) hDFs rapidly filled cell-free area in ~30 hours upon making a scratch, while Oct4+ hDFs were much slower filling cell-free area in ~50 hours (Figs 1e and S3), indicating that Oct4 over-expression impairs cell migration.
Altogether, the observed change of Oct4+ cells morphology, changes in the levels of mesenchymal-and epithelial-related markers, and slower cell migration might suggest that hDFs undergo mesenchymal-to-epithelial transition (MET) during Oct4-induced cell plastic state. miRNA-Seq results: sample to sample variation and quality check. MiRNA expression was analysed using three independent biological replicates represented by three different hDF cell lines expressing Oct4 or GFP respectively (here referred to as hDF1-3 Oct4 or hDF1-3 GFP).
At day 6 post transduction and antibiotic selection, total RNA was isolated from control GFP+ and Oct4+ hDFs (see Fig. 2a for the experimental design) and subjected to miRNA-Seq. Every biological replicate contained more than 4.5 × 10 6 non-filtered reads and Cook's distance analysis did not reveal any outliers among sequenced biological samples (Fig. S5). Hierarchical clustering, PCA analysis, and correlation matrix between samples showed highly distinct miRNA expression profiles between Oct4+ and GFP+ hDF cells, while there was no significant intra-group variation from sample to sample (Fig. 2b-d).
www.nature.com/scientificreports www.nature.com/scientificreports/ miRNA-Seq results: differentially expressed miRNAs. Given the striking difference in miRNA expression profile between Oct4+ and GFP+ hDFs, we sought to further characterize miRNAs in cell plasticity induction process. We detected 1,654 miRNAs from approximately 2,600 human mature miRNAs described in miRBase database 31 , while 212 were significantly (p < 0.05) differentially expressed between Oct4 and GFP expressing hDFs. When more stringent criteria (p < 0.05, log2(fold change) > 1.5) were applied, then 42 miR-NAs were significantly up-regulated and 20 down-regulated in Oct4+ hDF when compared to GFP controls (Fig. 3a,b).
The expression of several differentially expressed miRNAs (up-regulated: miR-96-5p, miR-182-5p, miR-183-5p, miR-9-5p, miR-302a,b,c,d-3p; down-regulated: miR-34c-5p, miR-193a-5p, miR-424-5p, miR-503-5p, www.nature.com/scientificreports www.nature.com/scientificreports/ miR-143-3p, miR-145-5p) was confirmed using RT-qPCR demonstrating that the expression of selected miRNAs was in concordance with the miRNA-Seq data (Fig. 4a,b). It is of note that the expression of some miRNAs is different among the replicates. This may be caused by different levels of Oct4 expression among the replicates, as the expression of some miRNAs (e.g. members of mir-302 cluster) correlates with the expression level of Oct4 in individual hDF cell lines (Figs 1b and 4a). Additionally, triplicates represent three different cell lines derived from three different individuals, thus the variability of miRNA expression might be also contributed to the nature of distinct hDF cell lines.
To address the question, whether changes to the miRNAs expression upon Oct4 overexpression is temporal or permanent, we assessed the expression of key selected miRNAs at day 6 and 20 upon transduction (Fig. S6). Although the expression of Oct4 gradually increases, reaching ~900 fold change at day 20, not all miRNAs shared www.nature.com/scientificreports www.nature.com/scientificreports/ similar expression pattern, indicating that Oct4-induced plastic state is not stable at miRNA expression level. Some miRNAs from mir-302 cluster (miR-302a-3p, miR-302b-3p) were significantly up-regulated at day 20, while miR-193a-5p and miR-145-5p were significantly down-regulated, when compared to day 6. miRNA-Seq results: genomic clusters and family analyses. It becomes increasingly clear that a single unique miRNA can be associated with specific cell functions. However, frequently miRNAs act in concert, being www.nature.com/scientificreports www.nature.com/scientificreports/ transcribed as clusters and/or acting as miRNA families, targeting the same target(s) or pathways(s). Therefore, better understanding of the function of the miRNAome is given rather by analysing the whole miRNA groups than individual miRNA molecules. We performed clustering of all expressed miRNAs based on their genomic cluster position (miRNA clusters, arbitrarily defined in miRBase as miRNAs located within 10,000 bases) and seed sequence similarities (miRNA families). We found 19 clusters and 15 miRNA families to be up-regulated and www.nature.com/scientificreports www.nature.com/scientificreports/ 10 clusters and 13 families to be down-regulated in hDFs over-expressing Oct4 when compared to GFP controls (p < 0.05, fold change > 1.5) (Fig. 5a-d).
In many scenarios, the genomic locations of differentially expressed miRNAs are not random, but they are clustered into specific loci of the genome. It is possible that these loci share the same regulatory pathways, however little is known about them. We visualised these loci into Circos plot, which shows a disperse distribution of the loci across the genome, with chromosome 1 being the most enriched with detected differentially expressed miRNA clusters followed by chromosomes 2 and 7 (Fig. 5e). Interestingly, many differentially expressed miRNA clusters are also miRNA families (red lines in the centre area), for example: mir-302, mir-182, mir-192, and mir-34.
As Oct4 over-expression in fibroblasts promotes cell fate switch into the neural and hematopoietic lineages, we sought to compare the list of up-regulated miRNAs with miRNAs expressed in neural and hematopoietic lineages. We used miRmine database 34 that provides a global view on human tissue miRNA expression profiles data. We found that 12.5% (18 miRNAs) of upregulated miRNAs were also found to be expressed in brain, 7.64% (11 miRNAs) are found to be expressed in blood, and 13.89% (20 miRNAs) were expressed in both brain and in blood (Fig. 6b).

Over-expression of mir-302 cluster in hDFs alters cell morphology and down-regulates mesenchymal markers N-cadherin and β-catenin.
Mir-302 cluster is the most up-regulated miRNA cluster in Oct4+ hDFs (Figs 3a and 5a) and it is classically viewed as a pluripotency-associated miRNA cluster that promotes reprogramming to iPSCs and it is also capable of altering cell fate 35 . We over-expressed mir-302 cluster in hDFs using Doxycycline-inducible vector (Fig. 7a). Similarly to Oct4 over-expression, cells changed their morphology upon mir-302 over-expression with transition of long-spindled fibroblast morphology to short-spindled cell shape that was associated with altered migration capacity of cells, and down-regulated mesenchymal markers N-cadherin and β-catenin (Fig. 7b-d).

Discussion
Reprogramming of somatic cells to iPSCs is of great research importance and high hopes have been pinned on the potential application of iPSCs in regenerative medicine. However, time requirements and high costs of the reprogramming process, as well as possible tumorigenicity, immunogenicity and increased genomic instability of the resulting cells represent the major limitations towards practical application of iPSCs 36,37 . Therefore, an alternative reprogramming approach involving cell plasticity induction has been proposed.
Cell fate conversion using Oct4-mediated cell plasticity induction is a complex process regulated by a large number of gene networks and signalling pathways 4,5 , potentially containing also miRNAs. MiRNAs, a family of small non-coding RNAs that regulate gene expression, represent the key regulators of virtually all processes in cells and it is well documented that manipulations with single or multiple miRNAs can influence altering cell fate (reviewed in 21 ). Still, the role of miRNAs in this process is not fully understood.
Here we show that Oct4 expression in hDFs down-regulates mesenchymal-related, up-regulates epithelial-related genes, and alters cell migration, suggesting that MET might play a role in this process. MET is crucial for early stage of the reprogramming process towards the generation of iPSCs 38,39 , however it is not clear whether MET is crucial for Oct4-mediated plasticity. In the first study describing the Oct4-induced plastic cells, MET has not been observed presumably due to different culture media, containing bFGF and IGF2, used in this study 5 . However, the authors observed profound changes in Oct4+ fibroblasts to acquire cuboidal shape that should represent a predictor for subsequent cell conversion competency, thus the involvement of epithelial/mesenchymal-associated genes in this process remains to be determined. It has been also demonstrated that Oct4-induced plasticity is characterized by mixed expression of lineage-specific genes 5 . We have not found any significant up-regulation of selected neuro-or cardio-associated genes (Fig. S4). However, we did not use media containing linage-inducing growth factors in our experimental setup, thus the mixed expression of lineage-specific genes may require the presence of these growth factors.
Our small RNA-Seq analysis revealed that Oct4+ hDFs possess distinct miRNA expression profile when compared to controls. We found more than 200 differentially expressed miRNAs upon Oct4 expression in hDFs. Interestingly, the data shows similarities with our previous results from miRNA microarray analysis comparing hDFs and iPSCs 25 . Among the top significantly down-regulated miRNAs clusters in hDFs over-expressing Oct4 are mir-424/503 and mir-34 clusters, whose members (e.g. mir-424, mir-503, mir-542, mir-34c), are found to be down-regulated in iPSCs as well 21,25 . We also found mir-143/145 cluster to be down-regulated upon Oct4 over-expression corroborating our previous findings that miR-145-5p is highly expressed in fibroblasts and inhibition of this miRNA in fibroblasts induces MET and facilitates reprogramming into iPSCs 25 . Surprisingly, we found embryonic-specific mir-302 and mir-183 clusters significantly up-regulated upon Oct4-induced plasticity in hDFs. Mir-302 cluster (comprised of mir-302a, mir-302b, mir-302c, mir-302d and mir-367) is referred to as embryonic stem cell-specific miRNA cluster that maintains pluripotency and ensures short G1 phase, therefore preserving unique properties of undifferentiated pluripotent cells 40 . It is also of note, that mir-302 cluster www.nature.com/scientificreports www.nature.com/scientificreports/ (alone or in combination with small molecules) has been shown to reprogram differentiated cells into iPSCs 41 . Interestingly, hDFs inducibly over-expressing mir-302 cluster show similar features as Oct4+ hDFs shortly upon induction -cuboidal shape, changes in the levels of mesenchymal markers and altered cell migration. www.nature.com/scientificreports www.nature.com/scientificreports/ One could argue that mir-302 cluster up-regulation in Oct4+ hDFs is not so surprising, as Oct4 directly activates its transcription in pluripotent stem cells 40,42 . However, in OSKM reprogramming mir-302 cluster is induced only late in the process, moreover in embryonic stem cells its expression requires Oct4 in orchestration with Sox2 in order to maintain pluripotent state 21,25,40 . Thus, mir-302 cluster up-regulation shortly upon Oct4 over-expression is rather unexpected, indicating a unique molecular profile of Oct4+ hDFs. Similarly, members of mir-183 cluster (mir-96, mir-182 and mir-183) are expressed at high levels in pluripotent stem cells and are gradually expressed only during the late stage of OSKM reprogramming 43,44 . On the contrary, we did not detect significantly elevated levels of other pluripotent stem cells-enriched miRNAs, such as mir-205, mir-17-92 cluster, mir-371-373 cluster 21,33 , suggesting that these miRNAs are not crucial for Oct4-mediated cell fate alterations and underscoring the unique miRNA expression profile of Oct4+ hDFs.
As we have detected some differentially expressed miRNA clusters to share the expression pattern with pluripotent stem cells, one could expect that Oct4-mediated cell fate switch requires a brief transition through pluripotent state. It has been demonstrated that brief exposure to OSKM factors and subsequent cell fate switch requires this transition 8,9 , however global gene expression profiling of Oct4+ fibroblasts did not reveal any up-regulation of genes associated with pluripotency or genes typically expressed during early or late stage of iPSCs reprogramming process 5 . This is in concordance with miRNA expression profile of Oct4+ hDFs indicating mixed miRNA expression, however biased towards embryonic-like pattern. Together these finding indicate that forced Oct4 expression in fibroblasts leads to unique cell plastic state associated with gene and miRNA expression profile that is distinct from cells undergoing reprogramming towards iPSCs.
For gene expression assays, the isolated RNA was reverse transcribed using Transcriptor First Strand cDNA Synthesis Kit with Anchored-oligo(dT) 18 Primer (Roche). Reverse transcription products were then amplified by quantitative PCR (LightCycler ® 480, Roche) using LightCycler ® 480 SYBR Green I Master. The primers used in this study are listed in the Fig. S7. Relative gene expression was determined using ΔΔC t method and normalized to endogenous control GAPDH.
NGS and library preparation. NEBNext Multiplex Small RNA Library Prep Set for Illumina (Set1, 2; New England Biolabs; Ipswich, MA, USA) was used to prepare libraries for further sequencing according to the manufacturer's instructions as described previously 49 . Briefly, 800 ng of total RNA were used to create size selected small RNA library. The sequencing (single end) was performed with 2.0 pM library using the NextSeq ® 500/550 High Output Kit v2 (75 cycles; Illumina; San Diego, CA, USA).
NGS data analysis. The data were processed as described previously 49 . Briefly, the quality of the raw sequencing data was assessed using FastQC (v0.11.5) 50 , Minion and Swan (Kraken package, v15-065) 51 were used to scan and identify adaptor sequences which were subsequently removed by Cutadapt (v1.12) 52 . Only adapter-containing reads were kept for further processing. The adapter-trimmed reads were further processed using the following steps: (1) Removal of very low-quality read ends (Phred < 5); (2) Keeping only reads with Phred score of 10 over at least 85% of the length; (3) Only reads within 16-28 bp were kept as potential microRNA reads. FASTX-Toolkit (v0.0.14) 53 was used for the quality filtering, the rest of the steps were performed by Cutadapt (v1.12) 52 and bash scripting. FastQ Screen (v0.10.0) 54 and Bowtie (v1.1.2) 55 were used to evaluate overall mapping of the preprocessed reads to the human reference genome (hg38) 56 and to determine possible contamination by rRNA and tRNA molecules. An additional quality check was performed by Sequence Imp (v15-075) 51 . The raw microRNA expression levels were quantified by Chimira (v1.0) 57 . All data analysis and visualizations were conducted using R 3.4.4 58 with installed Bioconductor project 59 . The following R packages were used to analyse and visualize data: DESeq2 60 , ggplot2 61 , pheatmap 62 , RColorBrewer 63 , plotly 64 , Mirbase.db 65 , and RCircos 66 packages were used to cluster and visualize miRNAs into families and clusters.
Scratch-wound assay. Cells were plated into 12-well plate (150,000 cells/well), cultured for additional three days and then scratched using a 100 µl pipette tip. Plates were photographed immediately after scratching (5 positions for each condition) and then every 1 hour for total of 50 hours using time-laps Cell R imaging station (Olympus). Cell-free area was measured using ImageJ (https://imagej.nih.gov/ij/) software.

Data availability
The datasets generated during and/or analysed during the current study are available in the GEO database repository (https://www.ncbi.nlm.nih.gov/geo/), accession number: GSE126284