Solexa Sequencing Identification of Conserved and Novel microRNAs in Backfat of Large White and Chinese Meishan Pigs

The domestic pig (Sus scrofa), an important species in animal production industry, is a right model for studying adipogenesis and fat deposition. In order to expand the repertoire of porcine miRNAs and further explore potential regulatory miRNAs which have influence on adipogenesis, high-throughput Solexa sequencing approach was adopted to identify miRNAs in backfat of Large White (lean type pig) and Meishan pigs (Chinese indigenous fatty pig). We identified 215 unique miRNAs comprising 75 known pre-miRNAs, of which 49 miRNA*s were first identified in our study, 73 miRNAs were overlapped in both libraries, and 140 were novelly predicted miRNAs, and 215 unique miRNAs were collectively corresponding to 235 independent genomic loci. Furthermore, we analyzed the sequence variations, seed edits and phylogenetic development of the miRNAs. 17 miRNAs were widely conserved from vertebrates to invertebrates, suggesting that these miRNAs may serve as potential evolutional biomarkers. 9 conserved miRNAs with significantly differential expressions were determined. The expression of miR-215, miR-135, miR-224 and miR-146b was higher in Large White pigs, opposite to the patterns shown by miR-1a, miR-133a, miR-122, miR-204 and miR-183. Almost all novel miRNAs could be considered pig-specific except ssc-miR-1343, miR-2320, miR-2326, miR-2411 and miR-2483 which had homologs in Bos taurus, among which ssc-miR-1343, miR-2320, miR-2411 and miR-2483 were validated in backfat tissue by stem-loop qPCR. Our results displayed a high level of concordance between the qPCR and Solexa sequencing method in 9 of 10 miRNAs comparisons except for miR-1a. Moreover, we found 2 miRNAs, miR-135 and miR-183, may exert impacts on porcine backfat development through WNT signaling pathway. In conclusion, our research develops porcine miRNAs and should be beneficial to study the adipogenesis and fat deposition of different pig breeds based on miRNAs.


Introduction
MicroRNAs (miRNAs), the small non-coding RNAs (typically 19-23 nucleotides), were firstly discovered in C. elegans [1] and play important roles in regulating post-transcriptional translation [2]. These small non-coding RNAs extensively exist in the cells of plants, animals and some viruses [3][4][5]. Increasing evidence demonstrates that miRNAs have evolved as a major class of generegulatory molecules critical for diverse biological processes, such as development, differentiation, proliferation and diabetes in eukaryotes [6][7][8][9].
miRNAs are initially transcribed from miRNA genes located mainly in the intergenic genomic region. By RNA polymerase II, they yield transcripts called primary miRNAs (pri-miRNAs), which are processed to generate miRNA precursors (pre-miRNAs) by a complex containing Drosha and then by Dicer to excise miRNA:miRNA-star (miR:miR*) duplexes. Of which one strand is more stable and preferentially incorporated into an RNA-induced silencing complex (RISC), while the other strand (miRNA*) is usually discarded [10]. The miRNA then guides the RISC to bind the UTRs and CDS of target mRNAs, where it downregulates the gene expression, often by blocking protein production or by degrading the target mRNA [11,12].
The pig (Sus scrofa), an important species in animal production industry, is a right model for studying adipogenesis and fat deposition. Increasing evidences suggest that miRNAs contribute to the regulation of fat development [13][14][15]. However, rare reports could be read about the porcine miRNAs that effect the adipocyte differentiation and adipogenesis. The total number of porcine miRNA genes discovered by human beings so far is much lower, which affects the study on pigs. Therefore, finding more miRNAs in pigs can be contribute to various and further studies on pigs. Some recent articles have depicted miRNAs in pigs [16][17][18][19], but most of these researches concentrated on miRNAs of porcine skeletal muscle. Backfat is another important tissue, however, there have been scarce reports on miRNAs identified from this tissue.
Solexa sequencing technology is a convincing strategy for identifying miRNAs and has been utilized by several laboratories to identify miRNAs [10,20,21]. To enrich the repertoire of miRNAs in pigs, we constructed and sequenced two small RNA libraries prepared from the backfat tissue of 150-day-old Large White (LW) and Meishan (MS) pigs. The former is a lean type pig, and the latter is a Chinese indigenous fatty pig. As the spatial distribution of miRNAs may contribute to the mechanistic understanding of adipocyte differentiation and adipogenesis, this study intended to provide deeper insights into dynamic regulation of the conserved and novel miRNA levels during different breeds of pigs.

Results and Discussion
Solexa sequencing data between the two small RNA libraries Small RNAs (18-30 nt) were obtained from the total RNA, 59 and 39 RNA adaptors were ligated to the RNA pool, and the adaptor-ligated small RNAs were then subsequently subjected to RT-PCR to produce sequencing libraries. PCR products were purified and small RNA libraries were sequenced using Solexa. Primary analysis was performed on the 35 nt tags. First, we got rid of the 59 adaptor, the low quality tags to obtain clean tags. The length distribution of the clean tags was subsequently summarized. Then, the clean tags were annotated into different categories, and finally those tags which couldn't be annotated to any category were used to predict the novel miRNAs. The work flow of Solexa sequencing is shown in Figure S1. The flow results of raw sequence data for the two libraries are shown in Table S1. The clean reads in the small RNA library of Large White was more than 19 million, and it was the same for Meishan small RNA library. Small RNA annotation is presented in Table S1, A-B. The distribution of sequence lengths was similar between both small RNA libraries by following the law of normal distribution, and the majority of sequences were 22 nt in length (Table S1, C-D). miRNA had a much tighter length distribution, centering on 21-24 nt. The number of 21-23 nt sequences was significantly greater than those shorter or longer sequences, and almost half of the sequences in Large White (49.79%) and Meishan (46.28%) are canonical 22 nt miRNA, which coincided with the known specificity for Dicer processing and the features of mature miRNAs [22,23]. Nearly 10 million reads of Large White and Meishan small RNA libraries were respectively screened as miRNA candidates and used for further analysis. The common and specific tags of two libraries, including the summary of unique tags and total tags were listed in Table S1, E in which the total number of unique sRNAs reads in Large White and Meishan small RNA libraries was almost 0.38 million and 0.49 million, respectively. These unique sRNAs sequences contained 0.1 million common sequences between the two libraries. The common sequences of total sRNAs between the two libraries were much higher (accounted for 97.75%) than that of unique sRNAs (accounted for 13.63%), this indicated that the number of Large White-specific and Meishan-specific sequences was much smaller than that of common sequences between the two libraries.
The known miRNAs expression profiles between two libraries are shown in Table 1. The expression of the known miRNAs in the two samples was demonstrated by plotting Log2-ratio and Scatter Plot (Figure 1). The results showed that 9 miRNAs were significantly different between the two libraries, such as the expression of miR-215, miR-135 and miR-146b was higher in Large White pigs, opposite to the patterns shown by miR-1a, miR-133a, miR-122, miR-204 and miR-183 ( Table 1), suggesting that these miRNAs may have effects on the development of backfat tissue. Detailed information of known miRNAs in both libraries, including sequence reads and hairpin structures are shown in Table S3. Lots of miRNAs were found to have end variants, differing by one or several nucleotides most commonly in the 39 end. For example, miR-204 and miR-450 only had 39 end variants. miR-107, miR-122 and miR-19a were respectively differed by only one nucleotide in the 59 end, but they had several 39 end variants. Such variants may come from altered miRNA processing, preferential degradation at the miRNA ends, or post-transcriptional modifications, including RNA editing [24,25].
miR-1a and miR-133a had been reported to be exclusively expressed in muscle [7]. In our analysis, we identified that they also expressed in porcine backfat tissue, the reads number of miR-1a and miR-133a was 0.2 million and 1028, respectively, indicating that these two miRNAs may be effective in adipogenesis. A number of miRNAs (miR-136, miR-140, miR-204, miR-221, and miR-325) originated from the 59 and 39 arms of the miRNA hairpin precursors, exhibited a similar number of sequence reads. The approximate equivalent expression rates of miRNA and miRNA* (miR-136:293, miR-136*:107 in Large White library) may be due to the similar 59 end stability which results to equal incorporation of either strand into the RISC and protection from degradation [26]. Accordingly, this miRNA precursor may produce functional molecules on both arms.

Identification of novel porcine miRNAs
Base on the Solexa sequencing, we identified 140 novel porcine miRNAs, which corresponding to 156 genomic loci. 87 in Large White pigs, 102 novel miRNAs in Meishan pigs. Their secondary structures are shown respectively in Table S4. 2 of the 156 genomic loci were located in the exons of protein-coding genes, 35 in introns, 11 in the introns and exons of single protein-coding genes, 3 in the introns and exons of two protein-coding genes, others were in intergenic regions (Table S2). Almost all these miRNAs could be considered pig-specific except ssc-miR-1343, miR-2320, miR-2326, miR-2411 and miR-2483 which had homologs in Bos taurus ( Table 2).
The read number for each novel miRNA was much lower than that for the majority of conserved miRNAs. For instance, ssc-miR-2411 and ssc-un0012, both were novel miRNAs, the total reads were only 131 and 5, respectively. However, the total reads of two conserved miRNAs (ssc-miR-148a and ssc-let-7c) were striking, almost 1.7 million and 1.6 million, respectively (Table S2). This might suggest a specific role for these novel miRNAs during developmental stages. Since this study investigated the miRNAs of backfat tissue from female Large White and Meishan pigs, whether these low-abundant miRNAs are expressed at higher levels in other tissues and organs remains to be investigated. Future experiments could provide more insight into the function of these miRNAs.

Validation of porcine miRNAs
We adopted stem-loop qPCR to validate our sequencing data [27]. A total of 4 novel miRNAs were cloned and validated in backfat tissue by qPCR, and they all could express (Figure 2, A). The expressions of miR-1343, miR-2320 and miR-2411 were slightly higher in Meishan samples; miR-2483 higher in large White samples (Table S2). And our qPCR results showed the consistency (Figure 2, B). In addition, the qPCR results suggested that 6 miRNAs which were sequenced by Solexa and displayed significantly differential expressions were expressed in porcine backfat tissue, but the expression levels of different miRNAs varied (Figure 3). In general, the results of qPCR validated the sequencing results and they were consistent. The expression of miR-135, miR-224 was significantly higher in Large White pigs, opposite to the expression of miR-133a, miR-122. The expression of miR-204, higher in Meishan pigs, was significant between the two samples ( Figure 3). However, miR-1a was an exception, its expression level was higher in Large White pigs, which was inconsistent with the sequencing results, and no significant difference was observed. The discrepancy in the miRNA expression may because of  differently adopted methods, such as cloning protocols, and various reads numbers that could be produced. Our results displayed a high level of concordance between the two methods in the 9 of 10 miRNAs comparisons, showing that the majority of our miRNA expression data embody the actual miRNA expression levels.

Seed edit of miRNAs and clusters of miRNA sequences
The position of 2,8 of a mature miRNA is called seed region which is highly conserved, and the seed region most often binds to a target site in the 39 UTR of the target mRNA by perfect complementarity. The target of a miRNA might be different with the change of nucleotides in this region and any nucleotide difference in seed region may lead to the dysfunction of miRNA and then resulted in the variation of phenotypic characteristics [28,29].
In our analysis, miRNAs which might have seed edit can be distinguished by matching unannotated sRNA with porcine mature miRNAs from miRBase15.0, allowing a mismatch on certain position. We found in the two samples a total of 22 mature miRNAs with single nucleotide substitution in seed region (Table  S5). Our study demonstrats that A -to-G, G -to-A, C -to-T, Gto-T, A -to-C and T -to-C were the dominant substitutions, which is consistent with published reports [30,31]. In our research, small RNA sequences not perfectly aligning to the genome were often discarded as mismatched sequences, which led to 'sequencing errors'. However, RNA editing or post-transcriptional modifications, not absolutely technical artifacts, could also be attributed to 'sequence errors' [30]. Obviously, high-abundant miRNAs (let-7c, let-7f, miR-148a, miR-21 and miR-24) had higher edited probability in backfat tissue. This indicates that highly-expressed miRNAs targeted more genes in this tissue. Accordingly, nucleotide changes in seed region could increase miRNA diversity, change their targets and may have important phenotypic consequences.
Of the conserved miRNAs, 66 miRNAs, almost occupy one-third of total amount, were also found in previous report [32]. However, the remaining miRNAs are distinct, this may be due to the different porcine breeds and developmental stages. Several miRNAs such as let-7, miR-27 and miR-103 are reported to perform very important functions in adipogenesis [33][34][35]. Moreover, published report showed that miR-196, miR-214, miR-199b, miR-186, miR-101 and miR-27a were related to bovine backfat thickness [36]. These miRNAs were also observed in our miRNAs expression patterns,  suggesting that they may exert impacts similar to those of their orthologs on porcine adipogenesis and adipose tissue.
The 215 miRNAs were categorized into clusters on the basis of their locations on chromosomes (miRNA-miRNA distance,10 kb). In brief, 18 miRNAs (13 conserved miRNAs and 5 predicted miRNAs) were found to be in 7 different clusters (Table S6). The mir-17-mir-20 cluster consisted of four miRNA genes. By using RNAfold [37], its folding structure was predicted (Figure 4), which was consistent with the previous report in human [38]. Such clustered miRNAs could be coordinately transcribed together as polycistronic transcripts.
We implemented the phylogenetic analysis of mir-17 cluster among different mammalian species using MEGA4 [39]. mir-17 and mir-20 had closer phylogenetic and evolutionary relationship with each other (Figure 5). The result partially proved that mir-20 resulted from duplication of mir-17 and the history of mir-17 cluster might be closely linked to the early evolution of the mammalian lineage [38]. Furthermore, we identified that 17 miRNAs were extensively conserved from vertebrates to invertebrates (Table S7), suggesting that these miRNAs may serve as potential evolutional biomarkers.

Potential regulation pathway of miRNA in porcine adipogenesis
In the present study, the fact that miR-122 expressed in backfat is inconsistent with current report that this miRNA is liver-specific [40]. Researches have shown that miR-122 regulates the pathway of cholesterol biosynthesis, promotes lipid synthesis and inhibits fatty-acid oxidation in the adult mouse liver. In the present study, several key genes known to regulate fatty-acid synthesis and oxidation were downregulated, including FASN and ACC after miR-122 inhibition and their fatty-acid synthesis rate was also reduced [14]. We therefore could postulate that miR-122 may influence lipid metabolism and adipogenesis in the abovementioned way in backfat tissue. However, the fine-tuned mechanism is worthy of further research.
Recently, a novel mechanism for APC regulation has been found. It was shown that miR-135 targets the 39 UTR of APC and suppresses its expression, then APC-free b-catenin stimulates the WNT signaling pathway, leading to active transcription of target genes, and induce downstream WNT pathway activity [41]. Previous report has demonstrated that WNT signaling pathway inhibits adipocyte differentiation by blocking the expression of C/ EBPa and PPARc, two transcription factors indispensable for adipogenesis [42]. So, we could speculate that miR-135 may suppress porcine adipogenesis through activating WNT signaling pathway by targeting APC (Figure 6).
To explore the potential function of miRNAs with significantly differential expression in backfat tissue, we adopted algorithms RNA22, RNAhybrid, TargetScan, DIANA LAB and PicTar to predict the targets of miRNAs. The predictions were conducted according to the interactions of human mRNA-miRNA because of the absence of porcine miRNAs in current versions of the abovementioned algorithms. We found LRP6, the surface co-receptor component of WNT signaling pathway, had putative binding site in its 39 UTR (Figure 7), and this indicated that LRP6 may be the potential target of miR-183. Previous research demonstrated that LRP6 co-receptor is critical for regulation of adipogenesis and potentially for development of white adipose tissue [43], suggesting that miR-183 may exert an impact on adipogenesis through targeting LRP6.
In conclusion, we have identified 75 known miRNAs and 140 novel miRNAs in backfat tissue from Large White and Meishan female pigs. The study expands the repertoire of porcine miRNAs and could be contributive to various and further studies on fat deposition and adipogenesis of pigs. In addition, it indicated that  several potential miRNAs may influence adipogenesis and fat deposition. However, whether these miRNAs we have identified also expressed in male pigs or other tissues, organs and whether the potential miRNAs really function in adipogenesis and fat deposition remain to be investigated.

Ethics statement
All study involving animals were conducted according to the regulation (No. 5 proclaim of the Standing Committee of Hubei People's Congress) approved by the Standing Committee of Hubei People's Congress, P. R. China. Sample collection was approved by the ethics committee of Huazhong Agricultural University. The approved permit number for this study is ''No. 30700571''. Animals were allowed access to feed and water ad libitum under same normal conditions and were humanely sacrificed as necessary to ameliorate suffering.

Animal and sample collection
Six purebred female pigs (3 Large White and 3 Meishan pigs, 150-day-old) provided by Jingpin Pig Station of National Engineering Research Center for Livestock (Huazhong Agricultural University) were used in the present study. The backfat tissues were harvested and immediately frozen in liquid nitrogen, and then kept at 280uC until subsequent analysis.

Preparation of small RNA library and sequencing
Backfat at the 11-12th rib of 150-day-old Large White and Meishan pigs were collected for RNA isolation. Total RNA was isolated using Trizol agent (Invitrogen), according to the manufacturer's instructions. A 1% agarose gel, stained by ethidium bromide, was run to preliminarily indicate the integrity of the RNA. All RNA samples were quantified and examined for protein contamination (A260 nm/A280 nm ratios) and reagent contamination (A260 nm/A230 nm ratios) by a Nanodrop ND 1000 spectrophotometer.
In this study, two miRNA libraries were constructed, total RNAs extracted respectively from three Large White and Meishan pigs were pooled. Total RNA were prepared for small RNA Sequencing-by-Synthesis according to the procedure and standards of the Illumina Sample Preparation Protocol.

Sequencing data analysis
The raw sequence reads were generated by the Illumina Genome Analyzer at BGI-Shenzhen, China.
In order to determine conserved miRNAs, the filtered sequences were initially used to search the miRBase 15.0 with BLASTN, which allowed a maximum of two mismatches, and the gaps were counted as mismatches. The criterion was then implemented according to reported miRNA protocol [44].
Potentially novel miRNAs were identified by folding the flanking genome sequence of unique small RNAs using MIREAP( https://sourceforge.net/projects/mireap/).
Comparison of the known miRNA expression between two samples was conducted to find out the differentially expressed miRNAs. The expression of miRNA was shown in two samples by plotting Log2-ratio figure and Scatter Plot. The procedures are shown as below: (1) Normalize the expression of miRNA in two samples (backfat of Large White and Meishan) to get the expression of transcript per million. When the normalized expression of a certain miRNA was zero between two samples, we revised its expression value to 0.01. While if the normalized expression of a certain miRNA was lower than 1, further differential expression analysis was conducted without this miRNA. Normalized expression (NE) = Actual miRNA count/ Total count of clean reads61000000. (2)  The x and y represent normalized expression level, and the N1 and N2 represent total count of clean reads of a given miRNA in small RNA library of backfat of Large White and Meishan pigs, respectively [32].

miRNA validation via stem-loop qPCR
Stem-loop qPCR method [27] was used to validate the conserved and novel miRNAs. For each miRNA, three biological replicates were performed. Briefly, the assay was performed using stem-loop RT-PCR followed by SYBR Green Real-time PCR analysis. The DDCt method was used to determine the expression level differences between samples [45]. The significant level was set to 0.05. For the RT-PCR, we adopted the RevertAid TM first strand cDNA synthesis kit (Fermentas, K1622), according to the manufacturer's instructions. Porcine U6 snRNA was used as an internal control and all reactions were run in triplicate. The mature miRNA and primer sequences are available in Table S8.

miRNA targets predicted by bioinformatics method
To explore the potential function of miRNAs with significantly differential expression in porcine backfat tissue, the targets of miR-183 were predicted by algorithms RNA22, RNAhybrid, TargetScan, DIANA LAB and PicTar [46][47][48][49][50]. Figure S1 The work flow of Solexa sequencing. A. The experiment process. B. The whole data analysis process. (TIF)