Ethics Statement

Texel and Ujumqin sheep were obtained from a sheep stud farm located in Youyu, Shan Xi Province. All experimental and surgical procedures were approved by the Biological Studies Animal Care and Use Committee, Shanxi Province, Peoples Republic of China. The ewes were housed in one group and were fed according to the nutrient requirements of sheep established by the National Research Council in 1985; the feeding was in line with the Instructive Notions with Respect to Caring for Laboratory Animals that was published in 2006 by the Science and Technology Department of China (Approval No. S20072911).

Sample Collection and RNA Extraction

Five Texel and eight Ujumqin ewes of similar age (3–5 years old), body weight (50–55 kg), and body size were selected. The estrus of the 13 ewes was synchronized and artificial insemination was completed. The date of artificial fertilization was used as day zero of gestation. Texel and Ujumqin fetuses were selected from five stages in utero by surgery: 70, 85, 100, 120, and 135 days after mating, abbreviated below as 70 d, 85 d, 100 d, 120 d, and 135 d, respectively. Ujumqin lambs were euthanized at three stages after birth: 0, 35, and 70 days, abbreviated below as 0 d, 35 d, and 70 d, respectively. The longissimus dorsi muscles (LM) of these 13 animals were dissected and quickly stored in liquid nitrogen.

Total RNA was separately isolated and purified from 13 frozen LM samples using a special Animal Tissue RNA Purification Kit (Product #TRK-1002; Qiagen, Hilden, Germany). After purification, RNA purity was assessed using a spectrometer (Nanodrop, Wilmington, DE, USA) at ratios of OD₂₆₀/OD₂₈₀ above 1.8. This showed that RNA was high purity. And RNA integrity assessing was conducted using agarose gel electrophoresis (AGE). Two 28S/18S rRNA bands were seen clearly in AGE without other bands. The results illustrated that RNA was high integrity and could be further used for Solexa deep sequencing.

Library Construction and Deep Sequencing

For Texel sample set, equal quantities of total RNA isolated from the longissimus dorsi muscles of Texel fetuses at 70, 85, 100, 120, and 135 d were pooled (see Table S1). For Ujumqin sample set, equal quantities of total RNA isolated from the longissimus dorsi muscles of Ujumqin fetuses at 70, 85, 100, 120, and 135 d, and lambs at 0 (birth), 35, and 70 d were pooled. Two small RNA libraries were generated from equal quantities of total RNA representing Texel and Ujumqin sample sets using Illumina Truseq Small RNA Preparation kits. The two purified cDNA libraries were used for cluster generation with Illumina’s Cluster Station and then sequenced on an Illumina GAIIx instrument. Raw sequencing reads were obtained using Illumina’s Sequencing Control Studio software version 2.8 (SCS v2.8) following real-time sequencing image analysis and base-calling using Illumina’s Real-Time Analysis version 1.8.70 (RTA v1.8.70). The extracted sequencing reads were used in data analysis.

Data Analysis

We used a proprietary pipeline script, ACGT101-miR v4.2, to analyze the sequencing data. The data analyses were performed as follows:

1. Obtaining Mappable Sequences from Raw Sequencing Data. “Impure” sequences due to sample preparation, sequencing chemistry and processes, and the optical digital resolution of the sequencer detector were removed using a series of digital filters. The remaining sequencing sequences (sequ seqs, between 15 and 32 nt in length) were grouped into families (unique sequences).
2. Mapping miRNA-mappable Unique Sequences (Unique Seqs) to Pre-miRNAs and the Genome. The number of selected mammalian pre-miRNAs in the latest version of miRBase (v17.0) was 6,752, among which the number of Ovis aries pre-miRNAs was 55. Unique seqs were aligned against pre-miRNAs of Mammalia.
3. Mapping Unique Sequences to Selected Databases.

For conciseness, the classification of six generated sequence groups was summarized in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Classification of six generated sequence groups by data analysis.

https://doi.org/10.1371/journal.pone.0058905.t001

Prediction of Ovis aries miRNA Target Genes

In principle of miRNA target with mRNA, we downloaded 3536 Ovis aries mRNA sequences from NCBI database. Then we extracted 2258 conserved 3′UTR sequences of Ovis aries genes. Applied target prediction tool which was TargetScan database (http://www.targetscan.org/) to predict the potential Ovis aries target genes of 89 differentially expressed miRNAs between Texel and Ujumqin sheep (see Table S2). A series of rules were used: (1) a perfect Watson-Crick match between target gene and miRNA at 2–8 positions (numbered from the 5′ end); (2) one G:U pair in the seed match was allowable; (3) the threshold for the minimum context score percentile of the seed match was 50.

Quantitative Real-time PCR of Ovis aries miRNAs

The pool RNA used in deep sequencing was reverse transcribed by stem-loop antisense primers (see Table S3), respectively. The expression of 25 selected miRNAs was determined using the Quantitative real-time PCR conducted on the ABI PRISM® 7900HT Real-time PCR Detection System, with 5S ribosomal RNA as the internal control gene. For cDNA synthesis, 250 ng total RNA and 0.5 µL RT primer (2 µmol/L) were mixed and ddH₂O was added to adjust the total volume to 3.0 µL, the mixture was incubated at 65°C for 10 min and snapped on ice for 3 min. Then a 5.0 µL reaction mixture contained 3.0 µL denatured total RNA and RT primer (2 µmol/L), 0.25 µL dNTP (10 mmol/L each), 1.0 µL 5×RT buffer, 0.25 µL RNase inhibitor (40 U/µL) and 0.5 µL M-MLV (200 U/µL) and reaction conditions were 42°C for 60 min, and 70°C for 15 min and hold at 4°C. The quantitative real-time PCR (qRT-PCR) reaction mixture (20 µL) contained 0.5 µL RT product, 0.8 µL primer mix (10 µmol/L), 10.0 µL 2 X SYBR Green Mix With ROX, 8.7 µL ddH₂O and cycling conditions were 50°C for 2 min; 95°C for 2 min; and followed by 39 cycles of 95°C for 15 s, and 60°C for 30 s. Melting curves were analyzed after amplification. All reactions were run in triplicate. We analyzed relative quantification results by the 2 ^−ΔΔCT method, and statistically significant differences between Texel and Ujumqin were examined using a t-test.

Solexa Sequencing Data

Raw sequencing reads were obtained using Illumina’s Sequencing Analysis software; 35,700,772 raw reads were generated in the two libraries. We used ACGT101-miR v4.2 to analyze the sequencing data (see Figure S1). By removing various un-mappable sequencing reads from the raw sequence reads, or the sequ seqs, the remaining sequ seqs (filtered sequ seqs with length between 15 and 32 nt) were further classified as mappable sequences (see Figure S2). After collating and mapping sequences to the constructed databases (sheep genome project and Mammalia miRBase), overlapping reads were gathered into clusters. Reads within each cluster were sorted by the number of mismatched bases, and counts were blasted to the genome. The same mismatched reads were arranged in line with the expression of the reads.

The 2,319 unique miRNAs were divided into three groups: high count miRNAs (counts ≥4,718), middle count miRNAs (10≤ counts <4,718), and low count miRNAs (0≤ counts <10). Each library included 2,319 miRNAs that represented 10,940,530 total average counts (see Table S4). High count miRNAs accounted for 97.99% of the total miRNA expression on average, but only for 4.87% of the types of miRNA, suggesting that a few miRNAs played a leading role as candidate miRNAs. The average counts for high count miRNAs reached 94,878.89, whereas the average counts for low count miRNAs were only 1.38, indicating that only a few types of miRNA comprised the majority of sequences.

In our data, RFams (RFam: rRNA, tRNA, snRNA, snoRNA, and others) represented the maximum proportion of the known classes of RNA sequences (3.23%) and accounted for 5.99% of the total variation. However, mRNA represented only 0.12% of the total number of sequences and only 2.05% of the total variation (see Figure S3 and Table S5). The proportion of long-chain RNAs was very low, demonstrating that our constructed RNA libraries were of high quality and that the reliability of the sequence results was high. Nohit group reads were un-mapped to selected pre-miRNAs (in miRBase v17.0), mRNA, Rfam, Repbase, or the sheep genome. Because the sheep genome was not complete to date, 30.27% reads were classified as Nohit group reads. The sheep genome project is currently in progress and will lay the foundation for miRNA mapping from Nohit group reads.

Mappable Sequences

After excluding all un-mappable sequencing reads, the remaining sequ seqs were termed mappable sequences. The analysis of mappable sequences illustrated that the length distribution peaked at 21 and 23 nt, and 22 nt sequences accounted for the maximum percentage (51.69%). These results are typical products of Dicer incisions (see Figure 1A) [12]. Similar variation in miRNA length has been repeatedly detected in deep sequencing results of other species. The average phred score per base in a read after a 3′ adapter (3ADT) cut at base 37 accounted for the largest proportion of sequences (60%) (see Figure 1B). Both of these results confirmed that the sequ seqs were of high quality.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. The length and count distribution of the sequenced sequences.

A. Length distribution of reads after 3ADT (3′adapter) cut. B. Histogram of the average phred score¹ per base in a read after 3ADT cut. ¹ Phred score larger than 30 stands for probability of incorrect base calls less than 1 in 1,000 (above 99.9% accuracy) in one sequencing read. C. Histogram of the top 21 high count miRNAs (X-axis) versus their cumulative % (Y-axis) in total average counts of the two small RNA libraries. The cumulative % of the top 21 high count miRNAs and the overall 2,319 unique miRNAs were shown as the dashed horizontal line at 85.06% and the dashed horizontal line at 100%, respectively, and the % of the individual miRNA was expressed as the black line.

https://doi.org/10.1371/journal.pone.0058905.g001

Classification of Ovis aries miRNAs

Variation in the mappable sequences was concentrated at both the 3′ and 5′ ends, which produced multiple mature variants that were named isomiRs, as previously described [13]. Measurements of the abundance of a miRNA/miRNA* using the sum of the total isomiR sequence counts correlated well with the expression level of the most abundant miRNA/miRNA* sequence [14]. Thus, in our study, we focused on the most abundant miRNA/miRNA* sequence for the differential abundance analysis.

First, sequences were mapped to the constructed databases (see Figure S1). Then, sequences were sorted by their mapping to known ovine miRNAs and other known mammalian miRNAs, and their locations within the ovine genome. For conciseness, the seven types of obtained Ovis aries miRNAs were summarized in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of the seven types of obtained Ovis aries miRNAs.

https://doi.org/10.1371/journal.pone.0058905.t002

In total, 2,319 pre-miRNAs, expressing 2,914 miRNAs, were detected in this study. Some mature miRNAs might have been from different pre-miRNAs or genome locations, and therefore, these miRNAs were suffixed behind the name (see Table S20) according to the naming rule in miRBase. The analysis of the seven types of sequenced miRNAs was summarized in Tables S5. Ninety-six known Ovis aries miRNAs accounted for 11.67% of the total counts; 1,685 miRNAs were mapped to the sheep genome but were not homologous with other mammals. We thus identified them as Ovis aries-specific miRNAs; these miRNAs only accounted for 0.43% within the 11.67% of the total counts. The remaining 1,246 miRNAs could be mapped to other mammalian pre-miRNAs, and thus were described as mammalian conserved miRNAs; the portion of these miRNAs in the total counts exceeded 52.35%. These results also indicated that miRNAs were highly conserved in this species.

Analyses of sequence results revealed that the total counts were dominated by a few miRNAs (see Figure 1C); e.g., the first nine miRNAs accounted for 75.84% of the total average counts. In addition, the first nine miRNAs were all mammalian conserved miRNAs; two of the nine were known Ovis aries miRNAs, while the other seven were other mammalian conserved miRNAs.

The expression of 2319 miRNAs between Texel and Ujumqin sheep was further log₂-transformed. After analysis, 89 Ovis aries miRNAs were of significant difference (| log₂ | ≥1.5) and high level (counts ≥10) in the expression (see Table S2). These 89 miRNAs were selected for further prediction of miRNA target genes.

Known Ovis aries miRNAs

In total, 4,079,965 sequences were detected in gp1a, accounting for 11.67% of the total counts. These mappable reads were mapped to the Ovis aries genome and 103 Ovis aries miRNAs in miRBase v17.0 (see Figure 2A), and they were further classified as 96 miRNAs (see Table S6). These results indicate that miRNAs that were known and had higher expression levels were much easier to detect.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. The Ovis aries known miRNAs.

Comparison of Ovis aries known pre-miRNAs between miRBase v17.0 and gp1a in this study. A. The overall 55 Ovis aries pre-miRNAs recorded in miRBase v17.0. B. The 49 Ovis aries known pre-miRNAs sequenced in gp1a of this study. These pre-miRNAs in four groups: type containing miR-5p (tawny), miR-3p (light green), both miR-5p & miR-3p (aqua), and six undetected pre-miRNAs (lilac) in gp1a.

https://doi.org/10.1371/journal.pone.0058905.g002

Forty-nine of 55 Ovis aries pre-miRNAs in miRBase v17.0 were sequenced in gp1a (see Figure 2B), and they coded 89 known miRNAs. The high coverage illustrated that the constructed RNA libraries contained the majority of known Ovis aries miRNAs. After further study, we discovered that 41 of 49 detected pre-miRNAs coded both the 5′ and 3′ ends of miRNAs (miR-5p and miR-3p). In many cases, miRNAs that have rapid turnover rates cannot be detected by conventional methods, although deep sequencing allows many of them to be identified.

The expression patterns of the 103 known Ovis aries miRNAs differed considerably in this study, ranging from a few to millions of copies. Oar-miR-127-3p and oar-miR-379-5p were the first two miRNAs and they accounted for 5.13% and 2.60% of the total mappable reads, respectively. Six new sheep miRNAs were discovered in gp1a: oar-miR-127-5p, oar-miR-136-3p, oar-miR-323b-5p, oar-miR-323c-5p, oar-miR-431-3p, and oar-miR-432-3p. They showed complementary base pairing with known Ovis aries miRNAs (see Table S6).

isomiRs

The mappable reads in each group are listed based on the number of mismatched bases with known Ovis aries miRNAs, and reads of the same mismatched base number are also listed by the number of copies. For example, 6,849 isomiRs were detected for both the 5′ and 3′ ends of 49 pre-miRNAs in gp1a; within those, oar-mir-379 had 388 variants at the 5′ end and 96 variants at the 3′-end, respectively (see Table S7). These results fully illustrate that known miRNAs were mainly present as isomiRs. Moreover, 45 isomiRs that had the maximum counts in 96 miRNA clusters were consistent with known Ovis aries miRNAs, and another 51 isomiRs differed by very few bases from known sheep miRNAs.

The number of isomiRs of individual miRNAs ranged from one to hundreds; e.g., oar-mir-127 had the highest diversity of isomiRs (716), and isomiRs in 24 of 49 clusters in gp1a had more than 100 variants. The variation in isomiRs was quite diverse; differences occurred at only the 5′ end, only the 3′ end or at both ends. Variation at the 5′ end might affect the seed sequence (the 2^nd to 7^th bases at the 5′ end) and could produce different miRNAs, thereby changing the target mRNA, or cause upregulation of the transcription in some regions [15].

Candidate Ovis aries miRNAs

The expression levels of miRNAs discovered in this study are summarized in Table S20. We found that 94 miRNAs that corresponded to 93 unique miRNAs were known Ovis aries miRNAs. Another 2,820 predicted Ovis aries miRNAs, representing 2,226 unique miRNAs, were candidate Ovis aries miRNAs. Within these, 1,653 predicted candidate miRNAs, representing 1,216 unique miRNAs, could not be mapped to any known miRNAs, but they were mapped to the sheep genome and occurred within hairpins. This “PC” type miRNAs might represent Ovis aries-specific miRNAs. The other 1,167 miRNAs, representing 1,010 unique miRNAs, were mammalian conserved miRNAs.

The let-7 family of miRNAs comprises one of the key regulatory elements in the developmental process. The let-7 family has been studied in many species, including mammals, birds, insects, and plants. Phylogenetic analyses have demonstrated that the let-7 family is highly conserved in both sequence and function among mammals, and it plays a critical role during animal development. For example, let-7 was identified as a heterochronic switch gene. The loss of let-7 gene activity causes reiteration of larval cell fates during the adult stage, whereas increased let-7 gene dosage leads to precocious expression of adult fates during larval stages.

To date, nine kinds of let-7 family genes have been identified [i.e., let-7a, let-7b, let-7c, let-7d, let-7e, let-7f, let-7g, let-7i, and let-7j (only identified in dogs)] in mammals. Eight kinds of Ovis aries let-7 family genes were discovered during this research: oar-let-7a-2 (427,369.75 reads), oar-let-7a-3 (288 reads), oar-let-7b (84,793.63 reads), oar-let-7c-1 (86,997.63 reads), oar-let-7d (28,816 reads), oar-let-7e (18,198.73 reads), oar-let-7f-1 (156,359.98 reads), oar-let-7g (19 reads), and oar-let-7i (148,896.89 reads). None of these miRNAs were included in miRBase v17.0. Obvious differences were noted in the expression levels of the let-7 family. A comparison of these eight let-7 family types between Ovis aries and 12 other mammals illustrated that the let-7 family sequences had high similarity within Mammalia (see Figure S4). This study also found that the Ovis aries let-7 family miRNAs possessed the same seed sequence (5′-GAGGTA-3′).

Chromosomal Mapping of Sheep miRNAs

The sheep whole genome sequencing project is ongoing and the size has been estimated at 2.97 Gb. The sheep genome database (February 2010) released the base sequences of 26 pairs of autosomes and the X chromosome. After we mapped the 2,319 Ovis aries pre-miRNAs found in this study to all of the sheep chromosomes (see Figure 3A), we found that 1,754 pre-miRNAs were located on chromosomes and 1,879 unique miRNAs that accounted for 81.03% of the expression were mapped to 2,436 genome locations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. The chromosomal map of Ovis aries miRNAs.

A. Chromosomal locations of Ovis aries miRNAs: miRNAs on “+ strand” (red line) and on “− strand” (blue line) across 27 pairs of chromosomes and chrUN. B. Histogram of pre-miRNAs on the X chromosome (white bar graph) and autosomes (gray bar graph) sequenced in Ovis aries and published in other six well-studies mammals. The chromosomal information of pre-miRNAs in these six mammals was from Sanger miRBase v18.0.

https://doi.org/10.1371/journal.pone.0058905.g003

The number of miRNA loci on each chromosome was unequal; oar-miR-493-3p was only detected at one genome location on chromosome 18, whereas PC-3p-1818407_1 had the highest number of genome locations (23) and occurred on multiple chromosomes. An analysis of the density distribution of miRNA loci on 27 pairs of Ovis aries chromosomes (number of miRNA loci per Mbp of individual chromosome) indicated that the density distribution ranged from 0.46 to 2.40 miRNA loci per Mbp. In all, 2,436 miRNA loci were mapped to the Ovis aries genome (1,201.95 Mbp in NCBI), producing an average density of 2.03 miRNA loci per Mbp. Chromosome 6 (129.07 Mbp, possessed 60 miRNA loci) had the lowest density distribution at 0.46 miRNA loci per Mbp. Chromosome 18 (72.49 Mbp with 174 miRNA loci) had the highest density distribution at 2.40 miRNA loci per Mbp. The longest chromosome l (299.84 Mbp), and the shortest chromosome 24 (45.32 Mbp), had 212 and 64 miRNA loci, respectively, producing respective densities of 0.71 and 1.41 miRNA loci per Mbp. The number of miRNA loci on a chromosome’s positive- and negative-strands was approximately equal, at 1,327 and 1,109, respectively.

X-linked miRNAs

Most mammals apply the XX/XY sex-determination system and no miRNAs are present on the Y chromosome. The Ovis aries X chromosome is 129.14 Mbp in size, ranking 4^th largest among all of the chromosomes. In this study, 133 pre-miRNAs were located on the X chromosome, accounting for 7.58% ( = 133/1,754) of the genome mapped pre-miRNAs (see Figure 3B). And 141 unique miRNAs were located on the X chromosome, representing 150 genome loci that accounted for 6.16% of the total miRNA loci. The density distribution of miRNA loci on the X chromosome reached 1.16 miRNA loci per Mbp, which ranked 4^th in the density distribution of miRNA loci among all of the chromosomes. The proportions of pre-miRNAs on the X chromosome within the total genomes of six other mammals (in Sanger miRBase v18.0) are the following: primates: human (hsa: Homo sapiens) 7.09%, chimpanzee (ptr: Pan troglodytes) 10.72%; model animals: pig (ssc: Sus scrofa) 14.58%, mouse (mmu: Mus musculus) 10.90%, rat (rno: Rattus norvegicus) 13.93%. However, in the cow (bta: Bos taurus), the proportion is only 5.61%. The distribution of miRNAs among chromosomes showed higher densities of miRNAs on X chromosomes compared to the average densities on autosomes in eight mammalian species [42], which was consistent with prior studies that demonstrated their resistance to meiotic sex chromosome inactivation [43].

Prediction and Analysis of Ovis aries miRNA Target Genes

8737 Ovis aries miRNA target genes were predicted (see Table S22), according to the series of rules which were described in “Materials and methods”. Seen from the predicted results, one Ovis aries miRNA could regulate one or more target genes, i.e., bta-miR-214 predicted 319 Ovis aries target genes; and one Ovis aries target gene could be regulated by one or more miRNAs, i.e., EF462423.1 had 36 target sites of miRNAs.

Quantitative Real-time PCR Validation of Ovis aries miRNAs Expression

To validate the Solexa sequencing results, 25 miRNAs (significant difference (| log₂ | ≥0.5) and high level (counts ≥50) in the expression between Texel and Ujumqin sheep, and equal distribution in the seven types of obtained Ovis aries miRNAs) were selected for further analysis. These 25 Ovis aries miRNAs were conducted stem-loop Quantitative real-time PCR assays from the same RNA preparations used for the Solexa sequencing. All of these 25 miRNAs could be identified in Texel and Ujumqin sheep. As shown in Table 3, eight out of ten miRNAs showed similar expression patterns as those revealed by our Solexa sequencing analysis. For unknown reasons, the expression levels of bta-miR-451 and PC-3p-14023_92 were inconsistent with the Solexa sequencing results. The results confirmed that these 25 miRNAs exist in Ovis aries.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Quantitative real-time PCR results of 10 miRNAs in Texel and Ujumqin sheep.

https://doi.org/10.1371/journal.pone.0058905.t003

High Coverage of Sequencing Results

With the application of high-throughput sequencing technology, miRNA arrays can be used to identify known miRNAs, and yet they cannot detect unknown miRNAs. Deep sequencing, which takes advantage of both clone sequencing and bioinformatic prediction, is the primary method for detecting miRNAs. Deep sequencing can both detect low expression miRNAs and predict unannotated miRNAs. The Solexa Genome Analyzer, which employs Sequencing by Synthesis (SBS), is the typical second-generation sequencing instrument.

MiRNA sequencing of two small RNA libraries was performed using a Solexa Genome Analyzer. The large number of sequenced sequences ensured not only the high quality of sequence calibration, but also the comprehensive coverage of miRNA expression. The two sequence libraries included 35,700,772 raw reads that represented 34,964,457 mappable reads (15–32 nt in length; see Figure S2). Yao et al. obtained a total of 1,147,787 high-quality reads. The length distribution peaked at 22 and 23 nt, which is consistent with the results commonly expected for miRNAs [16]. The length variation was mainly affected by enzyme modifications, including RNA editing [17], 3′ editing [18], and exonuclease activity [19]–[20]; 64.07% of the mappable sequences, which accounted for 0.25% of the variation in the sequences, were identified as miRNAs or miRNA candidates (see Table S5).

Ovis Aries miRNAome

MiRNAs are posttranscriptional regulators of gene expression; they are new targets for revealing the molecular mechanisms that form traits. Combining high-throughput sequencing and bioinformatics, this study examined 2,914 mature miRNAs representing 2,319 unique miRNAs. Among them, 1,103 unique miRNAs were conserved within Mammalia and 1,216 unique miRNAs were Ovis aries-specific. The count distribution of sheep miRNAs was very uneven; the first 21 miRNAs accounted for 85.06% of the total average counts (see Figure 1C). Muscle-specific miR-l, miR-133, and miR-206 all showed high expression (40,757≤ counts ≤2,699,554) in this study. Highly expressed miRNAs tended to be more stable, whereas low expression miRNAs were easily influenced by development stages. This phenomenon might relate to changes in cell regulation.

MiRNA genes are first transcribed to primary transcripts (pri-miRNAs) [21]–[25] that are further incised to about 70–90 nt hairpin precursors (pre-miRNAs). Then pre-miRNAs are processed into mature single-stranded miRNAs. Both the 5′ and 3′ ends of oar-mir-329a and oar-mir-329b matched miRNAs in miRBase v17.0, but only their 3′ ends (miR-3p) were detected in gp1a. This result agrees with previous reports in which at least 40% of miRNA sequences that were deposited in miRBase (v8.2) did not represent the predominantly cloned sequences, either because of differences in 3′ or 5′ end processing, or because of miRNA/miRNA* strand selection [26].

This study also found that the let-7 family of miRNAs possessed the same seed sequence (5′-GAGGTA-3′) reported in previous studies. Requiring conserved Watson-Crick pairing to the 5′ region of the miRNA centered on nucleotides 2–7, called the miRNA “seed”, markedly reduced the occurrence of false-positive predictions [27]–[30]. The 5′ region was the most conserved portion of metazoan miRNAs and the 5′ region of certain Drosophila miRNAs matched perfectly with 3′ UTR elements mediating mRNA decay and translational repression [31]. Subsequent experiments revealed that miRNA-like regulation was most sensitive to nucleotide substitutions that disrupted seed pairing [32]–[34].

Research on the miRNAome is the first step to understanding the expression, regulation, and function of miRNAs. Here, we produced a detailed overview of the Ovis aries miRNAome for the first time. Few papers have reported the identification of sheep miRNAs to date. These papers used bioinformatics, routine techniques, or deep sequencing based on partial tissues, and discovered 181 Ovis aries miRNAs [35]–[38]. Among them, 100 miRNAs (55.25%) were identical to our sequenced miRNAs and only seven miRNAs (3.87%) were not detected in our study. The remaining miRNAs had small differences from our sequenced miRNAs (see Table S21). This overlap in detected miRNAs among studies might be due to the high-expression levels of these miRNAs. However, other miRNAs showed high specificity: some miRNAs were only expressed during specific development stages, as part of certain diseases, or in specific species.

Clusters of miRNAs are also expressed as long primary transcripts. Located in a polycistron, co-expressed miRNA clusters are pivotal in coordinating the regulation of multiple processes, including embryonic development, the cell cycle, and cell differentiation [39]–[40]. By targeting components that have different roles along a signaling pathway, different members of one miRNA cluster can act as a whole to control the signal transduction process [41]. Although a large number of miRNA clusters has been discovered in animal and plant genomes, the functional consequences of this arrangement remain poorly understood.

isomiR

Further analysis of our sequence results revealed special characters of the Ovis aries miRNAome (e.g., almost all of the sequenced miRNAs produced different isomiRs), illustrating that the known miRNAs are mainly present in the form of isomiRs. In this study, for a given kind of miRNA, the corresponding number of isomiRs ranged from one to hundreds, and no direct relationship was observed between the expression of a miRNA and the kinds of isomiRs, which is similar to the findings of other deep sequencing studies. Kuchenbaue et al. reported that the number of isomiRs showed only a moderate correlation with the absolute expression levels of each miRNA (R² = 0.40, Pearson’s correlation coefficient), suggesting that the number of observed isomiRs was not directly related to the abundance of a miRNA. In some cases, the counts of isoforms were higher than the counts of corresponding known miRNAs in miRBase v17.0. We suggest that these most frequent isoforms should be used to refine miRBase annotations of Ovis aries miRNAs.

Every isomiR family has a feature and spatiotemporal specificity in certain physiological processes. More and more studies have revealed that isomiRs function in animals and may be marks of specific biogenesis processes and/or functions. This study has offered the most complete and accurate list of ovine isomiRs to date. The comprehensiveness and effectiveness of this information will set the stage for robust research into the complex functions of these regulatory molecules, which is necessary to decode the Ovis aries miRNAome in detail. The presence and various functions of isomiRs still require further research.

The English in this document has been checked by at least two professional editors, both native speakers of English. For a certificate, please see: http://www.textcheck.com/certificate/wi5a3Z.