Integrating miRNA and mRNA Expression Profiling Uncovers miRNAs Underlying Fat Deposition in Sheep

MicroRNAs (miRNAs) are endogenous, noncoding RNAs that regulate various biological processes including adipogenesis and fat metabolism. Here, we adopted a deep sequencing approach to determine the identity and abundance of miRNAs involved in fat deposition in adipose tissues from fat-tailed (Kazakhstan sheep, KS) and thin-tailed (Tibetan sheep, TS) sheep breeds. By comparing HiSeq data of these two breeds, 539 miRNAs were shared in both breeds, whereas 179 and 97 miRNAs were uniquely expressed in KS and TS, respectively. We also identified 35 miRNAs that are considered to be putative novel miRNAs. The integration of miRNA-mRNA analysis revealed that miRNA-associated targets were mainly involved in the gene ontology (GO) biological processes concerning cellular process and metabolic process, and miRNAs play critical roles in fat deposition through their ability to regulate fundamental pathways. These pathways included the MAPK signaling pathway, FoxO and Wnt signaling pathway, and focal adhesion. Taken together, our results define miRNA expression signatures that may contribute to fat deposition and lipid metabolism in sheep.


Introduction
Adipocytes are cells that form and store fat globules in the body. Adipose tissue deposits are located in different regions of the body, primarily in the form of subcutaneous fat or intramuscular fat in domestic animals. Though the volume of adipocytes varies from breed to breed and is influenced by diet, genetic factors are considered major determinants for the formation of adipocytes. Fat-tailed sheep, which exhibit distinctive large tails and hindquarters, comprise approximately 25% of the world's sheep population [1] and are raised commercially for meat, milk, fat, or wool production. Fat-tail is regarded as an adaptive response to the harsh challenges of desert life and the fat deposits provide a valuable energy reservation during drought seasons. The Kazakhstan sheep (KS, fat-tail breed) and Tibetan sheep (TS, thin-tail breed) are two native breeds raised in the extremely arid regions of western China, exhibiting distinct phenotypes of tails as previously described [2].
MicroRNAs (miRNAs) are a set of small noncoding RNAs (∼22 nt) that bind to partially complementary sequences on target mRNAs resulting in posttranscriptional regulation of gene expression [3]. MicroRNAs play key roles in regulating numerous biological processes in developmental, cell differentiation, and disease processes. These include adipogenesis and obesity in humans [4][5][6][7], as well as lipid metabolism and fat deposition in livestock [8,9]. Additionally, evidence suggests that miRNAs regulate the formation of adipose tissue. For instance, miR-378 was strongly associated with the thickness of back fat in beef cattle [10]. In recent years, next generation sequencing has been widely used to explore molecular mechanism of adipogenesis in sheep. However, miRNA-mRNA correlation analysis relating to adipogenesis in sheep has not been examined as in other livestock species [11,12].

BioMed Research International
This research provides an understanding of the process in sheep and contributes to the understanding of obesity and related diseases.
We have previously described specific changes of transcriptomic differences with respect to adipose tissues from the representative KS and TS breeds and identified a list of candidate genes underlying the phenotypic differences of fat/ thin tails in sheep [2]. In this study, we conducted miRNA sequencing using the same tissues, in attempt to elucidate the roles of miRNAs in the presentation of these two extreme phenotypes. We further employed a systems biology approach to examine the correlation between miRNA and mRNA expression data to identify potential miRNA-associated target genes. Here, we have implicated interactions from particular cellular processes such as the MAPK signaling pathway, FoxO and Wnt signaling pathway, and focal adhesion, gaining insights into the potential roles of miRNAs in the regulation of these pivotal pathways.

Animals and Phenotypes.
Six adult individuals (three males and three females, two years old) from KS and TS breeds were randomly selected, respectively, and any two or more individuals with a traceable phylogenetic relationship were avoided in the sampling process. Adipose tissue biopsies were collected from the tails and stored as previously described [2].

Small RNA Library Construction and Sequencing.
Total RNA from the mixed adipose tissues of six adult sheep was isolated using the RNAiso plus kit (TaKaRa, Dalian, China) according to the manufacturer's protocol. The small RNA libraries were prepared following Illumina5 TruSeq6 Small RNA Sample Preparation protocol. The small RNA libraries were sequenced using an Illumina/Solexa 1 G Genome Analyzer System (BGI, Shenzhen). The sequencing reads were deposited to the Sequence Read Archive, under accession code SRP093866.

Expression
Profiling. The raw miRNA sequencing reads were filtered and only the unique mapping reads were used for subsequent bioinformatics analysis, such as annotation and gene expression. The expression of miRNA in two samples (TS and KS) was normalized by transcripts per million (TPM) as previous description [13]. After normalization, the expression value of a miRNA, more than one in the two samples, was kept for the differential expression analysis, and we revised 0 value to be 0.01 when we calculated fold change of the two breeds [14]. The differentially expressed genes (DEGs) were determined according to its absolute fold change (log 2 (TS/KS)) of ≧1 and value of ≦0.05. Scatter plots were used to demonstrate differentially expressed miRNA between the two sheep phenotypes.

Quantitative RT-PCR.
Total RNA from the mixed sheep tails (six from KS and six from TS) was isolated using a miRNA extraction kit (CWBIO, Beijing, China). First strand cDNA was synthesised with reference to RevertAid First Strand cDNA Synthesis kit (Thermo Scientific #K1622, Fermentas) manufacturer's instructions, using special stem-loop RT primers [15] (Additional file 1 in Supplementary Material available online at https://doi.org/10.1155/2017/1857580). PCR-iQ5 system (Bio-Rad, Hercules, CA, USA) was used to analyse miRNA genes expression level. The reaction mixture volume is 25 L, containing 2 L of cDNA (25 ng), 12.5 L SYBR Premix Ex TaqTM II (TaKaRa, Dalian, China), 1.0 L specific forward primer, 1 L universal primer, and 8.5 L water, the reaction program referenced to previous description [14]. U6 snRNA was used as an endogenous control. All reactions were performed in triplicate, and 2 −ΔΔCt method was used to analyse the data. -tests (parametric test, unpaired test) was used to determine significance of the results using SPSS21. Graphpad Prism 6 was used to plot the graphs. All data were presented as mean ± SD and values of < 0.05 were considered statistically significant.
2.6. mRNA Profiling Bioinformatics. The sequences of reads with poly(A) tracts were deposited to the Sequence Read Archive, under accession code SRP093860. The FastQC package was used to assess the quality of raw reads and then mapped to the sheep genome assembly Oar v3.1 (GCA 000298735.1) using TopHat2 [16]. Quality control was performed on the two samples. The high-quality (HQ) clean reads were presented after a series of procedures (remove adapter, ploy N, and low-quality reads). All the following analyses were based on the HQ clean reads. The unannotated HQ clean reads were assembled by Cufflinks [17]. Orthonormal expression levels were calculated as FPKM (Fragments Per Kilobase of transcript per Million mapped reads) [18], and FPKM distribution of the samples was presented as gene's coverage. The edgeR package for bioconductor was applied to identify differentially expressed genes between TS and KS [19]. Genes were defined as differentially expressed with the threshold FDR < 0.05 and |log 2 (fold change)| > 1.

Integration Analysis of miRNA and mRNA Sequencing Data.
MicroRNA-mRNA integration analysis was performed on the mRNAs that were differentially expressed between the KS and TS breeds in our previous study [2], and the miRNAs were identified and classified according to the miRNA family again, in order to meet the current study. Individual miRNAs determined to be differentially expressed were submitted to each of the three databases ((RNAhybridv2.1.2)+svm light(v6.01), Miranda(v3.3a), and TargetScan(Version: 7.0)) to determine the predicted mRNA targets, and the parameter was taken as the default. We selected only mRNA targets that were differentially expressed in the same samples. Next, the correlation between the expression levels of each miRNA and any of the predicted, differentially expressed target mRNA from any of the three databases was computed. According to the analysis results and the prior research, we chose miRNA-125a-5p and one of the target genes (ESRR ) for validation analysis.

Functional Enrichment of Target Genes.
For each miRNA with a significant association with its set of mRNA targets (based on the gene set enrichment test), the GO and KEGG terms were extracted for all of its mRNA targets. We reveal the functions significantly related to predicted target gene candidates of miRNAs used gene ontology (GO) analysis. The primary pathways of the candidate target genes were determined by KEGG pathway analysis. The formula used to calculate this is the same as that used in the GO analysis. Genes with < 0.05 are considered as significantly enriched in target gene candidates. The KEGG indicates the main pathways in which the candidate genes are involved.

Small RNA Diversity in Sheep Adipose
Tissues. The sequencing of two small RNA libraries from the KS and TS yielded 12 M counts of sequencing reads, respectively ( Table 1). All of the clean reads were aligned with the sheep genome (Oar v3.1) sequence using SOAP software. The length distribution and the distribution of miRNAs on the genome were analysed. The high-quality reads in both groups exhibited the canonical size range distribution that is common to mammalian miRNAs [3] (Figure 1(a)). The vast majority of the reads were approximately 21∼23 nucleotides (nt) in length, and ∼22 nt reads accounted for 50% and 39% for each of high-quality reads from the KS and TS breeds, respectively. This is in agreement with previous founding in other adipose tissues [20]. These results demonstrated the high reliability of using the small RNA-sequencing approach to obtain miRNA reads for further studies.
There were 8,000,311 (67.53%) and 6,720,288 (56.68%) reads matched to the sheep genome sequence in KS and TS libraries, respectively. All clean reads were classified and annotated by tag2annotation software (developed by BGI), aligned against the Rfam database (Ver. 10.1) and the miRBase (Ver. 19.0). There were 11,599,356 and 11,603,988 total conserved miRNAs reads, and 0.77% and 0.81% of clean reads were identified as potential novel miRNAs from KS and TS libraries, respectively, and more information was given in the additional file (Additional file 2).

Differential miRNA Expression in Sheep Adipose Tissues.
To identify conserved miRNAs in sheep, all sRNA sequences were mapped to known miRNAs in the miRBase 19.0 database. We presented a total of 815 miRNAs found in abundance in adipose tissues of sheep (Additional file 3). Of 815 identified miRNAs, 539 miRNAs were expressed in both breeds, whereas 179 and 97 miRNAs were specifically expressed in KS and TS, respectively. The top ten abundant miRNAs from KS, TS, and both groups are listed in Table 2. Of the shared miRNAs, 93 were upregulated and 33 were downregulated in the adult adipose tissues of sheep ( Figure 1(b)).
We next analysed the differentially expressed miRNAs in the two breeds. One hundred and seventy-five miRNAs exhibited similar expression levels (nonsignificant), 33 miR-NAs showed significantly reduced expression ( < 0.05) and 276 miRNAs showed an increase expression ( < 0.05) in KS, compared to TS (Figure 1(b)). MicroRNAs with similar expression patterns in different sample pairs were clustered together. Clustering analysis was based on the sample difference model by using Cluster software, and the results were viewed with Java Treeview. All differentially expressed miRNAs clustered together after six rounds of clustering ( Figure 2).
Novel miRNAs could be predicted through the characteristic hairpin structure of miRNA precursors. We predicted novel miRNAs by exploring secondary structure, the Dicer cleavage site, and the minimum free energy of the unannotated small RNA reads, which could be mapped to sheep genome sequences by using Mireap software. In total, 186,755 unannotated sequences were used to predict novel miRNAs. Thirty-five potential novel miRNAs were identified; 7 miRNAs were expressed in the two samples, whereas 18 and 10 miRNAs were specifically expressed in KS and TS, respectively (Additional file 4).

Validation of miRNA by qPCR.
In order to validate the reliability of miRNA expression data, quantitative RT-PCR was conducted to assess the expression of six upregulated miRNAs (miR-27c, miR-183, miR-224, miR-339a,  20 21 22 23 24 25 26 27  18 Length (nt)  Figure 1: (a) Length distribution of the clean reads based on total abundance and distinct sequences. (b) Differential expression of miRNAs between the adipose tissues from KS and TS breed. Each point represents a single miRNA. The and axes show the expression level of miRNAs in these two samples, respectively. Red points represent miRNAs with a ratio > 2, blue points represent miRNAs with 1/2 < ratio ≦ 2, and green points represent miRNAs with ratio ≦ 1/2; ratio = normalized expression in treatment/normalized expression in control. miR-2070-3p, and miR-101) and four downregulated ones (miR-27d, miR-106a, miR-204, and miR-301) among the differentially expressed miRNAs, using adipose tissues from both KS and TS groups. U6 was used as an internal control. Stable miRNA expression levels in both KS and TS tissues correlated well between qPCR and miRNA-seq data ( Figure 3).

GO and Pathway Enrichment Analyses of Differentially Expressed Genes.
To determine the function of differentially expressed genes, all of the DEGs in this study were mapped to terms in the GO database (http://www.geneontology.org/) using Blast2GO [21]. In order to assess the general functional characteristics of genes activity in sheep adipose tissues, we performed functional enrichment analysis of GO using DAVID software and pathway enrichment using specific KEGG terms [22]. A total of 1522 genes were categorized into the three main categories of GO classification. Figure 5 shows the GO classification of DEGs. KEGG pathway annotation showed that 4899 genes and 537 DEGs genes were annotated for 209 biological functions (Additional file 9). Figure 6 showed the top 20 of pathway enrichment. We found that several significantly overrepresented categories of GO biological processes were associated with fat metabolism and deposition. We found that the cAMP signaling pathway had a critical role in both adipogenesis and lipid partitioning in white adipose tissue and wellcharacterized mechanisms controlling adipocyte differentiation [23,24].

Identification of Potential mRNA Targets of miRNAs.
MicroRNA recognizes its target mRNA through binding to a seed sequence, which localizes on the 3 untranslated regions (3 UTR) of target mRNAs [25]. In an attempt to determine potential genes whose mRNAs might be targeted by particular miRNAs, we conducted an integrated analysis between the expression of specific miRNA and the expression of all the predicted miRNA target genes (mRNA) using TargetScan.
We identified potential miRNA target genes that are downregulated at the transcriptional level and are inversely correlated with the miRNA expression in the same KS and TS adipose tissues. Table 3 lists the top 10 differentially expressed miRNAs and their target genes determined by miRNA-mRNA correlation analysis, and the details were presented in Additional file 10.

Validation of the Target Gene of miR-125a-5p.
For validation, we selected the a priori gene ESRR , which inhibits preadipocyte differentiation [26,27]. Using RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid) and miRanda (http://www.microrna.org) and combining the results of the analyses, we presumed that miR-125a-5p regulates ESRR translation. MiR-125a-5p was found to directly target the ESRR 3 UTR sequence (Figure 7(a)). To verify this finding, the 3 UTR containing miR-125a-5p targeted site was cloned and inserted into the psiCHECK6-2 reporter plasmid (psiCHECK6-2-ESRR -3 UTR WT) and also the mutated 3 UTR (psiCHECK6-2-ESRRa-3 UTR-mutant). The two types were cotransfected with the miR-125a-5p mimics or negative control sample into 293T cells, respectively, and luciferase activity was measured after 48 h later. The luciferase activity of the miR-125a-5p group was significantly lower than that of the NC group ( < 0.01), whereas the mutated ESRRa 3 UTR exhibited increased luciferase expression (Figure 7(c)). Thus, the miR-125a-5p was confirmed as one of the ESRR translation targeted regulation sequences and the primer sequences used for this are listed in Additional file 1.

Functional Enrichment Based on miRNA-Associated Target Genes.
In order to assess the general functional characteristics of the miRNA-related activity in sheep adipose tissues, we selected the putative miRNA target genes derived from miRNA-mRNA correlation analysis and performed functional enrichment analysis of gene ontology (GO) using DAVID software and pathway enrichment using specific KEGG terms [22]. We found that several significantly overrepresented categories of GO biological processes were associated with biological functions, including environmental information processing, cellular processes, metabolism, and genetic information processing (Additional file 11). KEGG pathway annotation showed that 4,687 target genes were annotated for 226 biological functions (Additional file 12).
Most of these genes were involved in metabolic processes, regulation of biological process, and signal transduction. Several of these KEGG terms were related to fat deposition and fatty acid metabolism, including MAPK signaling pathway, focal adhesion, Pyruvate metabolism, FoxO signaling pathway, and TNF signaling pathway. This suggests that

Discussion
Next generation sequencing approaches have become powerful tools to predict and identify the novel genes and small RNAs of livestock species and other organisms. Compared with other meat-producing livestock such as cattle and pig, fewer studies have focused on the role of miRNAs in the process of fat deposition in sheep, though there are a few omics studies screening miRNAs in sheep skin [28] and muscles [29]. In an effort to gain insights into the molecular mechanisms underlying fat deposition of adipose tissue governed by miRNA regulation in a fat-tailed and a thin-tailed sheep breeds, we evaluated whole transcriptome profiles of miRNAs and mRNA expression on the same adipose tissue from sheep tails [2].
In the present study, we were able to discover 35 novel miRNAs which have eluded previous efforts and identify 239 miRNAs that have not been previously annotated in sheep. Our analysis detected 815 differentially expressed miRNAs between the fat deposits of the two sheep breeds used here.
Not surprisingly, a proportion of the miRNAs was found to be differentially expressed in our study had been previously reported to play a role in adiposity, adipocyte development/differentiation, and other metabolic disturbances in livestock, as well as adiposity and related disease risk in humans. This further highlights the importance of miRNAs in adipose tissue development and metabolism. Previous studies have shown an increase in miR-143 during human and murine preadipocyte differentiation, and the following experiment verified that it inhibited preadipocyte differentiation [30][31][32]. miR-378 may promote bovine adipogenesis in white adipose tissue through targeting MAPK1 and PPAR [10]. The interactions demonstrated effects on the MAPK signaling pathway, focal adhesion, and Pyruvate metabolism pathways. miR-103 and miR-30 have been reported to enhance adipogenesis [32][33][34][35], while miR-27 and miR-138 could inhibit adipogenesis [36,37]. miR-122, miR-370, and miR-378 have been demonstrated to have a key role in lipid metabolism [7,38,39]. miR-148a might modulate fat deposition by targeting MAPKAPK5, MAPK3, and MAP2K2, and miR-148 has been shown to affect obesity and modulates adipocyte differentiation [40,41]. miR-125b-5p inhibited proliferation and promoted adipogenic differentiation in 3T3-L1 preadipocytes AGTGTCCAATTTCCC---------------AGAGTCCCT   Figure 7: miR-125a-5p downregulates the expression of ESRR by targeting its 3 UTR site. (a) The pairing schematic of ESRR 3 UTR and miR-125a-5p. The nucleotide in blue represents "seed sequence" of miR-125a-5p, and the mutation nucleotides are presented in red. (b) The insertion site of ESRR 3 UTR or its mutation luciferase reporter vector map. (c) The two types were cotransfected with miR-125a-5p mimics (or negative control) into 293T cells. Luciferase assay was performed 48 h after transfection. Results are presented as relative luciferase activity (mean ± SD; = 4; ANOVA; * * < 0.01). [42], and other studies have shown that miR-125b impaired brite adipocyte formation and function [43]. These reports are consistent with our findings reported.
In addition to microRNA expression profiling, we performed miRNA-mRNA integration analysis using the miRNA and mRNA expression data from corresponding tissues in distinct sheep breeds and determined a limited number of putative target genes for the differentially expressed miRNAs. The integration analysis removed a large portion of false positive miRNA target genes and increased the precision of our predictions [44]. Therefore, the identified candidates could be taken as reference genes to validate the differentially expressed miRNAs really affecting the mRNA expression levels during adiposeness. Here, we showed that miR-125a inhibits porcine preadipocytes differentiation by targeting ESRR [26]. This kind of targeted effect also exists in sheep, so we speculate that the function of miR-125a in ovine preadipocytes is similar to that in porcine preadipocytes.
The DAVID gene annotation analysis of the miRNAassociated targets indicated that the GO biological processes were enriched to modification processes, metabolic processes, signal transduction, and regulation of biological processes. This demonstrates the distinct roles of adipocyte development in the tail areas as the reserves of energy, in response to external or internal environmental conditions during the dry and cold seasons. KEGG pathway annotation of their putative miRNA-related targets revealed that the MAPK signaling pathway, focal adhesion, Pyruvate metabolism, FoxO signaling pathway, and TNF signaling pathway contribute to the formation of sheep fat/thin-tails phenotypes. Various tail phenotypes are unique for sheep compared to other animals in the word, interestingly, a mass of adipose depot in the fat-tail sheep breed and less in the thin-tail sheep breed. In a manner, this phenomenon is a form of energy storage in the fight against extreme conditions, and the function of tail adipose tissue is similar to subcutaneous adipose tissue of other animals according to the enrich pathways, which is associated with metabolism risk [20,45].

Conclusion
In conclusion, we identified a number of miRNAs that are differentially expressed between the fat-tailed and short-tailed sheep breeds. We further highlighted gene targets of related miRNAs that may be involved in regulating fat deposition and adiposeness, in sheep and other livestock. This occurs via the key signaling pathways including focal adhesion, Pyruvate metabolism, and the MAPK, FoxO, and TNF signaling pathway. Further studies are needed to verify the correlation between key miRNAs and their target genes by in vitro approach and elucidate the functional impacts that miRNAs serve during adiposeness. Our results also provide evidence for the interaction of miRNAs and genes in the regulation of obesity and metabolic syndromes, which suggests that this may serve as an animal model for human' obesity and metabolic syndromes researches.