Global transcriptome analysis of pig induced pluripotent stem cells derived from six and four reprogramming factors

Pigs are important, both for agriculture and as animal models for human diseases. However, due to the lack of embryonic stem cells, the possibility of genetic modification is quite limited. To overcome this limitation, induced pluripotent stem (iPS) cells have been derived from pigs. Despite the public availability of a large number of expression datasets from mice, rats, and primates-derived iPS cells, the expression profile of pig-derived iPS cells is quite limited. Furthermore, there is no dataset focused on the profiling of pig-derived iPS cell with six reprogramming factors (Oct3/4, Sox2, Klf4, c-Myc, Lin28, and Nanog). Here, we used Illumina RNA sequencing platform to characterize the mRNA expression of four-factor derived and six-factor derived pig iPS cells. We observed that the expression levels of whole genes in our established six factors derived iPS cells and parent fibroblast, and compared with that of iPS cells with four factors in public database. These data are valuable in understanding species difference in the reprogramming process of stem cells, and could help identify the key regulating genes involved in the process.


Background & Summary
The progress of next-generation sequencing technology has caused a technological breakthrough at the whole-genome level in a large number of species 1 . Especially, RNA-sequencing (RNA-Seq) has enabled us to take a snapshot of global gene expression in various organs and cells, regardless of any information of a reference genome. RNA-Seq outputs are digital data that can be uploaded to the public database, and sequence information can be shared worldwide.
RNA-Seq analysis also allows us to compare the biological similarity of embryonic stem cells (ES cell) with induced pluripotent stem (iPS) cells. In general, stem cells can be classified into two major subtypes: naïve and primed states 2,3 . ES/iPS cells at naïve state of pluripotency, reflect the characteristics of preimplantation embryos and are applicable in rodents, which contribute to chimeras and germ line 4,5 . The growth of naïve stem cell depends on the activation of LIF (Leukemia Inhibitory Factor) signaling, whereby the cell forms colonies with three-dimensional shape.
On the other hand, primed cells have the characteristics of post-implantation embryos and rarely contribute to chimeras and germ line. In brief, primed cells are already at a more differentiated stage compared to naïve cells. Primate ES/iPS cells were conventionally believed to be established in primed state 6,7 . However, recent publications have demonstrated a reliable method for transforming human ES cells from primed to naïve state 8,9 . Transcriptome analysis using RNA-Seq played an important role in identifying the cellular characteristics reported in those articles.
In the case of pig iPS cells, the status of the cells-naïve or primed-remains inconclusive since the pluripotent genes have a wide variety of phenotypes. To understand the biological variety of pig iPS cells, multiple datasets of global gene expression profiling would be needed. Although a significant number of reports on the establishment of pig iPS cells have been published [10][11][12][13][14][15][16][17][18][19][20] , expression profiling data in Sequence Read Archive (SRA) database are quite limited. Therefore, detailed biological features of pig iPS cells need to be addressed with whole expression profiling.
In our previous publication, we had reported that pig iPS cell, derived from six reprogramming factors, has more advantageous than that derived from four factors. Especially, the expression of six reprogramming factors was suitable for X chromosome re-activation 21 , which is one of the mile-stone characteristics of naïve-type stem cells. Our previous data using Ion Torrent sequencing also proved that the expression of six reprogramming factors was more advantageous to activate various pluripotent genes. Although the data obtained from Ion Torrent is suggestive, at least 20 M reads would be necessary to obtain a quantitative evaluation of the relatively low-expressing genes. The data obtained in our previous publication seem insufficient in terms of the number of sequencing reads required to conclude. This situation led us to detect the global expression profile of pig iPS cells, derived from the expression of six reprogramming factors, using Illumina short-read sequencer, HiSeq 2500. Currently, there are no publicly available dataset of six factor-derived pig iPS cells using Illumina sequencing platform.
The aim of this study was to clarify the difference of mRNA expression profiles between pig iPS cells derived from six and four reprogramming factors. We found relevant submitted data from two research groups on pig iPS cells with four reprogramming factors, in SRA 22,23 . We could compare ours with these gene expressions since both datasets were obtained with Illumina sequencer. In this study, we describe the detailed expression profile of pig iPS cells derived from four and six reprogramming factors. Multiple analyses demonstrated that the pig iPS cells derived from six factors formed independent clusters compared to those derived from four factors, and were distant from fibroblasts. Furthermore, we detected that the expression levels of various naïve-specific genes were relatively elevated in pig iPS cells derived from six factors. Our data set would contribute to the understanding of biological differences between the iPS cells derived from six and four reprogramming factors, and provide the scientific explanation of how diversity of pluripotency-related genes related to the process of animal evolution.

RNA preparation and sequencing
Pig-derived iPS cells were cultured on the feeder cells derived from mouse embryonic fibroblasts. When iPS cells reached to the confluent condition around 3 × 10 6 cells/35 mm diameter cell culture dish, the cells were lysed in 700 μl of the RA1 solution of NucleoSpin RNA extraction kit (code: 740955.250, Takara Bio, Shiga, Japan). The total RNAs from pig embryonic fibroblast, iPS clone 1 and iPS clone 2, three RNA samples were submitted to NODAI Genome Research Center for quality control and Next Generation Sequencing. All samples had integrity number (RIN) >9 in Bioanalyzer (Agilent, Santa Clara, CA, USA) 31 . TruSeq RNA libraries were prepared from total RNA according to the manufacturer's protocol (TruSeq RNA Library Prep Kit v2, Illumina, San Diego, CA, USA). Libraries for the three cDNA samples were sequenced using the Illumina Hiseq 2500 sequencing platform yielding about 45 million 100-bp paired-end sequence reads per sample.

Quality check and mapping
Quality of raw sequencing reads was evaluated using FastQC program. FastQC results plotted the position of nucleotide base call on the x-axis and their corresponding quality score on the y-axis. Figures  1 and 2 shows that all FASTQ sequencing files have mean quality score >36, which conforms with the manufacturer's standard >33 of our own data. The summary of the data were listed in Table 1. The sample information, such as single read, or paired-end, number of reprogramming factors for pig iPS,   were listed in Table 1. We analyzed the 9 data from reference and public data base, and 6 data from our own experiments. The raw read data were recorded in Data Bank of Japan (DDBJ), and have been assigned BioProject accession PRJDB5113 (Data Citation 1). Figure 3a shows our experimental workflow following quality validation of the sequence reads. All sequence reads were trimmed using PRINSEQ software for discarding the low-quality reads. In case the adapter sequence remained, cutadapt software was used. The sequencing reads were aligned to the Sus Scrofa Ensembl genome assembly (Sscrofa11.1) using STAR 24 . The rates of uniquely mapped reads ranged from 60-95% (Fig. 2b). After the mapping, the analysis was divided into two pipelines, cufflinks-R package and featureCount-TCC (Fig. 3a).

Counting gene expression and downstream analysis
For the cufflinks-R package pipeline, the distance matrix was obtained with cufflinks after the alignments. The differentially expressed genes were identified with R package, such as ape (dendrogram analysis), prcomp() in R (PCA analysis), Heatmap.3() in R (heatmap analysis) from results of cufflinks. For the featureCount-TCC analysis, the output with BAM format was processed with featureCount program, and the differentially expressed genes were identified with TCC-GUI program, which provided from website of Dr. Koji Kadota (Tokyo University, Tokyo, Japan) (https://infinityloop.shinyapps.io/TCC-GUI/). The basic algorism of TCC-GUI is TCC (an acronym for Tag Count Comparison) 25 and normalized gene expression with Edge R 26 . The parameter of TCC-GUI is described in below, normalization method was TMM method, differentially expressed genes were counted with EdgeR, the filtering threshold for low count genes was set at 7, FDR cut off was set as 0.1, elimination of potential differentially expressed genes was set at 0.05. The original code of analysis (two dimensional PCA plot, three dimensional PCA plot, dendrogram analysis, heatmap analysis, expression analysis of identified genes of TCC) were supplied as the function of TCC-GUI (https://infinityloop.shinyapps.io/TCC-GUI/). The output file for the normalized expression level with TCC-GUI (tmm_edger_3_0.1_0.05_TCC_Normalized.csv, Data Citation 2) in Figshare. The list of differentially expressed genes determined with TCC-GUI were also provided in Figshare (tmm_edger_3_0.1_0.05_TCC.csv, Data Citation 2).

Code availability
The following software and versions were used for quality control, data trimming, and data analysis as described in the main text:      Raw FASTQ file quality The results of FastQC showed that the quality score, as the accuracy of base call, was very high (mean > 36) compared to the Illumina standard >33 (Figs 1 and 2). Our sequences were mapped to the Sscrofa 11.1 genome assembly between 60-95% rates (Fig. 3b). After the mapping, analysis pipeline was divided into two lines, cufflinks-R package, and featureCount-TCC analysis. Results of two dimensional PCA plot were shown in Fig 3c and d. Although the two-dimensional PCA plots with both pipelines could not show the difference between fibroblast and pig iPS cells (Fig. 3d), three-dimensional PCA plot with featureCount-TCC pipeline showed the clear unique cluster of fibroblasts, pig iPS cells (Fig. 3e). Furthermore, the dendrogram analysis showed that cufflinks-R package pipeline divided into two groups, XX chromosome samples and XY chromosome samples (Fig. 4a). Interestingly, featureCount-TCC pipeline showed two clusters, fibroblast and pig iPS cells (Fig. 4b).
Comparison with published studies The current study was compared with two datasets from previous studies. We showed the expression of naïve-specific genes and primed specific genes as a heat map in Fig. 4c,d. As good agreement with dendrogram analysis, the heatmap analysis with featureCount-TCC analysis showed the unique cluster of pig iPS cells with six factor expression, while cufflinks-R package analysis could not form the unique cluster. We focused on three representative genes: ERAS, ESRRB, FDF4 which are well known as naïve markers. These genes had higher expression in our pig iPS cells with six factors than in pig iPS cells with four factors (Fig. 5a,b). These results indicate that pig iPS cells with six factors are located closer to naïve state.

Usage Notes
We have used Linux OS, Ubuntu 14.04.5 LTS for all analysis.