MicroRNA profiling of rhesus macaque embryonic stem cells

Isolation and characterization of rESC lines

Thirteen expanded rhesus macaque blastocysts with a prominent ICM were selected by immunosurgery and 11 ICMs were isolated and plated onto feeder cells. ICMs attached to feeder cells within 48 h and three ESC-like ICM outgrowths appeared after 7-8 days. ICM outgrowths were manually dissociated into 4-6 smaller clumps using a microscalpel, excised from feeder cells and replated onto fresh mouse embryonic fibroblasts (mEFs). Clones with distinct boundaries and high nuclear to cytoplasm ratios were selected for further propagation. Three rESC lines were established and designated as IVF1.2, IVF3.2 and IVF3.3. IVF3.2 and IVF3.3 were derived using the same sperm and oocyte donors. rESCs shared common morphologies with other primate ESCs such as being flat with a distinct boundary against feeder cells. Cells showed high nuclear to cytoplasm ratios and prominent nucleoli (Figure 1A). IVF1.2 and IVF3.3 were cultured for >60 passages and IVF3.2 for >80. Pluripotency markers were highly expressed in all rESC lines including Oct-4, Nanog, SSEA-4, TRA-1-60 and TRA-1-81 (Figure 1A). Other crucial transcription factors, such as Sox-2 and Rex-1, were also detected by RT-PCR (data not shown).

All rESC lines showed strong alkaline phosphatase activity (Figure 1A) and could spontaneously differentiate into cell lineages from the three embryonic germ layers. After 2 days in suspension culture, rESC clones formed embryoid bodies (EBs). After 5 days, cavities were observed in EBs, which are termed cystic EBs. EBs attached to gelatin-coated dishes and were cultured for 5-7 days until various differentiated cells appeared including neuron-like cells, contractive cardiomyocytes and endoderm-like cells (Figure 1B). Teratomas were formed 6-7 weeks after rESC injection into the hind leg muscles of SCID mice. Histological characterization revealed that ectoderm (neural tube and squamous epithelium), mesoderm (cartilage, muscle and adipose) and endoderm (intestinal epithelia and glands) like structures were present in teratomas (Figure 1B). These results demonstrated that rESC lines could differentiate into all three embryonic germ layers in vitro and in vivo.

Detailed G-banding analysis revealed that all rESC lines were karyotypically normal with a diploid set of 42 chromosomes (Fig. S1 in Additional file 1), even after long-term culture in vitro and repeated freeze-thaw procedures. IVF1.2 and IVF3.3 were derived from male blastocysts (40+XY) and IVF3.2 was derived from a female blastocyst (40+XX).

Small RNA Sequencing

Newly established rESCs were tested for purity using Oct-4 staining, which resulted in >95% Oct-4-positive cells (Fig. S2 in Additional file 1). A small RNA library from each cell line was then constructed for sequencing. Small RNA library sequencing yielded 12.66 × 10⁶, 13.12 × 10⁶and 11.57 × 10⁶ raw reads from IVF1.2, IVF3.2 and IVF3.3, respectively. After filtering, 10.89 × 10⁶ (IVF1.2), 10.60 × 10⁶ (IVF3.2) and 9.26 × 10⁶ (IVF3.3) clean reads (18-30 nt) were obtained. The length distribution of clean reads centered at 22-23 nt, coinciding with the miRNA length range. Sequence alignment with the rhesus macaque genome (rheMac2) using the Short oligonucleotide alignment program (SOAP) [25] demonstrated that 8.36 × 10⁶ (IVF1.2), 7.68 × 10⁶ (IVF3.2) and 7.15 × 10⁶ (IVF3.3) sequences were mapped to the genome (Figure 2A). MiRNAs were the major component of small RNA libraries from rESCs (Figure 2B) being 60% of sequences, which were matched to reference miRNAs in each sample. The proportions of other annotated RNAs was <15%. Notably, ~25% of sequences could not be annotated because of low coverage (6 × ) and relatively poor annotations of the available rhesus macaque genome. Unknown sequences were further analyzed to identify novel miRNAs.

MiRNA expression profiles of rESC lines

After matching to 466 reference rhesus macaque miRNAs, 326, 329 and 326 known miRNAs were detected in IVF1.2, IVF3.2 and IVF3.3, respectively, with absolute counts ranging from 1-1,702,656 (Table S1 in Additional file 1). Among them, 304 miRNAs were shared by the three pools and remaining miRNAs were expressed at low levels with <10 counts. From a genome-wide perspective, very low miRNAs were likely sequencing artifacts or miRNA machinery byproducts [26]. Therefore, an absolute count of 30 was set for each pool as the minimum threshold for further analyses.

To compare miRNA expression patterns between rESC lines, we evaluated the expression variability (Var) and coefficient of variation (C.V.) of miRNA clusters (Figure 3A, 3B). The results revealed similar expression patterns of miRNA clusters in rESC lines. We also compared the expression of total miRNAs using a differential index (D.I.) and Kappa Statistical analyses. For the D.I., most miRNAs were consistently expressed (Figure 3C), whereas a slight deviation existed in IVF1.2 due to low-level expressed miRNAs. For Kappa Statistical analysis, there was high repeatability between any two cell lines (Figure 3D). Taken together, miRNA expression patterns in rESC lines were highly consistent.

Identification of novel miRNAs

MiRNA target gene prediction and functional classification

Three programs were used to predict miRNA target genes. Using the TargetScan program, a custom set of Perl codes were downloaded from the TargetScan database (Release 5.1, http://www.targetscan.org/). Conserved and non-conserved predicted patterns of miRNA targets were considered separately. For PITA and miRanda, we used 238 expressed miRNAs (number of reads ≥30) and rhesus macaque 3'UTRs orthologous to annotated human 3'UTRs as inputs to predict miRNA targets. Over-representation of predicted miRNA targets for both conserved and non-conserved patterns based on the TargetScan program in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway is shown in Table S3 in Additional file 1 (P≤0.05, after Benjamini Hochberg correction). Gene enrichment patterns obtained by various methods were similar to each other, indicating that miRNAs correlated to functionally critical roles in rESCs. For example, predicted targets involved in the MAPK (hsa04010), Wnt (hsa04310) and TGF-beta signaling pathways (hsa04350), and pathways involved in cancer (hsa05200) were significantly enriched (Table S3 in Additional file 1).

MiRNA-target interaction network analysis

A systematic survey of miRNAs that interacted with targets may provide a better understanding of miRNA regulation in stem cells. Thus, we analyzed potential interactions between miRNAs and their targets, which is termed as the miRNA-target network (MT network). All miRNAs and their targets were used to generate a bipartite graph of miRNA-target interactions.. Nodes represent miRNAs or target genes and edges correspond to interactions between miRNAs with target genes. Node degree in the MT network is the number of connections or edges with other nodes. For each target gene node, degrees represent the number of miRNAs targeting the gene. Globally, the MT network was comprised of 8,934 nodes with 48,546 edges. Using MT network degree analysis (see methods), the degree of genes regulated by miRNAs was negatively correlated (r = -0.564 with p = 0.018 based on Pearson's correlation) with their evolutionary rate (dN/dS) (Figure 4A). Similarly, a previous study reported a negative correlation between the degree and evolutionary rate of proteins in a protein-protein interaction network [27].

To exclude factors with potential effects on the observed negative correlation, we tested the effect of 3'UTR length. After length normalization, the negative correlation remained (Figure 4B,4C). Another factor was a potential ascertainment bias caused by miRNA target prediction and identification of orthologous macaque 3'UTRs. To evaluate this possibly, we randomized the scale-free MT network based on Barabasi-Albert model. The degree of miRNA and target genes in the random network was significantly lower (p < 2.2e^-16 with a Wilcoxon rank sum test) compared with that observed in the MT network (Figure 4D). A significant correlation was not observed in the random network, which further supported the negative correlation in the MT network.

Comparison of miRNA profiles from mouse, human and rhesus macaque

To compare ESC miRNA profiles between species, we retrieved previously reported mouse and human ESC miRNA profiles. To avoid a cross-talking bias due to various sequencing techniques and data sizes, we selected Solexa sequencing datasets from mouse and human ESCs [14, 28]. MiRNA expression frequencies were normalized to each other by determining the expected frequency in mapped reads per million. MiRNA expression levels were calculated as the mean of three samples. The three species expressed 138 common miRNAs. Clustering analysis based on the 138 shared miRNAs revealed that miRNA expression profiles were more similar between human and rhesus macaque (Fig. S3 in Additional file 1), which was consistent with their evolutionary relationship. A number of miRNAs were stably and highly expressed between species, such as the miR-17 family, miR-21 and miR-103, which suggested conserved functions in ESCs during evolution.

Of the 138 common miRNAs the 'stemness' miR-290-295 cluster [13, 29] was the highest expressed in mouse ESCs. However, primate homologs the miR-290-295 and miR-371-373 clusters were not highly expressed in human and rESCs. Instead, the miR-302 cluster was the highest expressed in human and rhesus macaque [14, 30, 31], indicating the functional divergence of stemness miRNA clusters in primate lineages. The miR-290-295 cluster contains multiple mature miRNAs with seed sequences similar or identical to those in the miR-302 cluster [6]. MiR-290-295 and miR-302 clusters are transcriptionally regulated by Oct-4/Sox2 [29]. Taken together, these observations suggest that these two clusters are similarly expressed, regulate common target genes and are functionally conserved in ESCs derived from various species.

There are three miRNA clusters differentially expressed in mouse, human and rESCs. The chromosome 19 miRNA cluster (C19MC), a primate-specific miRNA cluster [32], was enriched in human ESCs, but almost absent in rESCs (Figure 5C). A conserved miRNA cluster within the imprinted Dlk1-Dio3 region was highly expressed in mouse and rhesus macaque ESCs, but rarely expressed in human ESCs (Figure 5B), which was consistent with results from four other human ESC lines [31].The miR-467 cluster in the sfmbt2 gene intron region was only expressed in mouse ESCs (Figure 5A), embryos and newborn ovaries [33, 34], but not in human and rhesus macaque ESCs, consistent with previous reports [29].

rESC isolation and characterization

rESCs were derived from the ICM of rhesus macaque blastocyst-stage embryos and characterized using standard protocols. rESC isolation and characterization details are provided in the supplemental Methods in Additional file 1.

RNA extraction and small RNA sequencing

Total RNA was extracted using a mirVana miRNA isolation kit (Ambion). According to the manufacturer's instructions (Illumina), the small RNA fraction between 18-30 nt was isolated from total RNA by PAGE (polyarclymide gel electrophoresis) purification and ligated to a pair of adaptors at the 5' and 3' ends. Small RNA molecules were converted to cDNA and amplified by RT-PCR using adaptor primers. Unique 6 bp indexed sequences were added to 5' adapters for multiplex sequencing. Purified DNA was directly used for cluster generation and sequencing analysis using an Illumina Genome Analyzer. These data have been deposited in the GEO database (GSE27886).

Small RNA annotation and novel miRNA prediction

For clean reads, the original 35 nt reads from sequencing were filtered by trimming the 5' and 3' adaptors, eliminating contaminants, and inadequate (<18 nt) and low quality reads. For detecting known miRNAs, clean reads were matched to reference miRNA precursor sequences in the rhesus macaque (miRBase version 16.0) using Blastall (-p blastn -F F -e 0.01). MiRNA absolute expression was the sum of clean reads mapped to the precursor. For annotating degradation tags of rRNA, tRNA,small cytoplasmic RNA (scRNA), small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA), clean reads were aligned with Genbank and Rfam 9.1 databases using Blastall. For identifying repeat associated RNA and degradation tags of mRNA, clean reads were overlapped with repeat sequences, exons and introns of mRNAs. In the initial alignment and annotation, some small RNA tags were mapped to more than one category. For uniquely mapped reads, the following priority rules were applied: rRNA/tRNA/scRNA/snRNA/snoRNA (Genbank > Rfam) > known miRNA > repeat > exon > intron. Unannotated small RNA tags were termed as 'unann'. Mireap (http://sourceforge.net/projects/mireap/) was used to predict novel miRNAs.

MiRNA expression profiling analysis

Based on the randomly selected 10,000 bins for each chromosome with each bin spanning a 50 bp window (210,000 bins in 21 chromosomes), miRNAs with sporadic reads were filtered out and miRNAs with >30 reads (p < 0.0001 based on Poisson distribution simulation) were treated as true reads for further analyses. miRNA expression profiles were compared and evaluated in rECSs at cluster and individual miRNA levels. A miRNA cluster was defined as consecutive miRNAs with a ≤10 kb inter-miRNA distance (MID). For each level, two independent methods were used for confirmation. Var and C.V. were used to survey miRNA member divergences in a cluster. Var was used to represent a divergence measure based on Shannon entropy [42]. D.I. and Kappa Statistical analysis (κ) were used to compare total expressed miRNAs in rESCs. D.I. was used to evaluate relative deviation of individual miRNA expression levels in rESCs. Kappa statistical analysis was used to measure the repeatability of total miRNAs expression between any two samples (supplemental methods in Additional file 1).

MiRNA target prediction

To predict miRNA target genes, TargetScanS [43], miRanda (http://www.microrna.org/microrna/home.do) and PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_prediction.html) were used. Common miRNAs (≥30 reads) were used to predict target genes. For each miRNA, the reference mature miRNA sequence in miRBase was used. For 3'UTR identification in the rhesus macaque, human annotated 3'UTRs from Refseq and orthologous 3'UTRs from six species (human, chimpanzee, macaque, mouse, rat and dog) were extracted from multi17way aligned files downloaded from the UCSC database [44]. Orthologous 3'UTRs were aligned using the MUSCLE [45] program to obtain multi6way 3'UTR files required for analysis.

Gene ontology (GO) analyses of predicted miRNA targets were performed using t DAVID Bioinformatics Resources [46] (http://david.abcc.ncifcrf.gov/home.jsp). Ascertainment biases such as a relatively low number of annotated miRNAs, 3'UTRs in the rhesus macaque, ≥30 cutoffs for miRNA reads and false predicted miRNA targets were assessed by simulating the bias with 240 randomly synthesized in silico 22 nt miRNAs (equal number with true expression).

Potential correlations between miRNA (or target) degrees connected in the MT network, miRNA expression and target evolution were evaluated (supplemental methods in Additional file 1).