Identication of the Genetic Basis for the Large-tailed Phenotypic Trait in Han Sheep Through Integrated mRNA and miRNA Analysis of Tail Fat Tissue Samples

Background: While evolution has led certain breeds of sheep to exhibit large tails composed of fatty tissue, the genetic basis for this fatty large-tailed phenotypic trait remains to be dened in breeds of Han sheep. Here, we employed a high-throughput sequencing approach to identify mRNAs and microRNAs (miRNAs) that were differentially expressed in tail fat tissue samples from large-tailed Han (LTH) and small-tailed Han (STH) sheep in order to identify key genetic determinants of the large-tailed phenotype. Results: In total, we identied 521 mRNAs (237 upregulated, 284 downregulated) and 14 miRNAs (6 upregulated, 8 downregulated) that were differentially expressed between these two sheep breeds. Predictive analytical database tools were subsequently utilized to identify 2,409 putative targets of these differentially expressed miRNAs (DEMs), including 65 which were among the list of differentially expressed genes (DEGs) identied in the present study. By specically focusing on predicted DEM/DEG pairs with appropriate regulatory directionality, we identied DIRF, HSD17B12, LPL, APOBR, INSIGI, THRSP, ACSL5, FAAH, ACSS2, APOA1, ACLY, and ACSM3 through mRNA analyses and ACSL4, FTO, FGF8, IGF2, GNPDA2, LIPG, PRKAA2, ELOVL7, SOAT2, and SIRT1 through miRNA analyses as candidate genes which may regulate fat deposition and fatty acid metabolism in the adipose tissue from the tails of Han sheep. Conclusion: Together, our data provide insight into the potential genetic basis for the large-tailed phenotype of LTH sheep, suggesting that it may be attributable to specic DEMs and DEGs that regulate one another and thereby control lipid metabolism. These data provide a basis for future research regarding the role of these genes in ovine tail fat deposition, and offer preliminary perspectives on the molecular mechanisms governing the fatty large-tailed phenotype in LTH sheep.

Multiple genomic regions liked to tail fat deposition have been identi ed in breeds of Iranian sheep [5].
Moioli et al [6] have also employed a genome-wide analytical approach which enabled them to identify genomic regions and functional genes in uencing ovine tail fat deposition. We have also utilized Ovine HD SNP chips to identify 8 quantitative trait loci associated with difference in tail fat deposition in Largetailed Han (LTH) and Small-tailed Han (STH) breeds of sheep through a range of statistical analyses (unpublished). Sample-size-weighted xed-effects genome-wide association (GWA) meta-analyses have further been used to clarify the genetic factors which in uence body fat deposition [4]. Recently, RNA sequencing has been employed as a means of assessing genome-wide transcriptomic activity and can be used to complement GWA analyses of complex traits [7]. In addition, microRNAs (miRNAs) can posttranscriptionally suppress mRNA expression, thereby further regulating phenotypes of interest [8]. Several studies have employed RNA sequencing strategies to clarify functional genes and gene networks regulating complex traits in humans [9], cattle [3], pigs [4,[10][11][12], and mice [13].
In an effort to identify genes linked to ovine tail fat deposition, Wang et al. (2014) analyzed the adipose tissue transcriptome of Kazak and Tibetan sheep breeds [14], while  conducted genome-wide analyses of mRNA and miRNA expression in STH and Dorset breeds to elucidate the genetic basis for phenotypic differences between these breeds [15,16]. While informative, these prior sequencing studies have focused on a limited number of animals and have not assessed intrinsic withingroup variability. Additionally, no prior studies have focused on integrative mRNA and miRNA expression pro ling in the adipose tail tissue of LTH and STH sheep breeds. As such, in the present study, we employed a high-throughput sequencing approach to identify mRNA and miRNA that were differentially expressed in the tail adipose tissue of LTH and STH animals with extremely divergent tail phenotypes. We further identi ed speci c differentially expressed genes and miRNAs (DEGs and DEMs respectively) linked to lipid metabolism, and performed an integrative analysis to elucidate the genetics basis for this large-tailed phenotype.

Animals sampling
In total, three unrelated adult LTH sheep and three unrelated adult STH sheep were selected from among domestic populations located in Jia county of Henan Province, China, and Cao county of Shandong Province, China. All experimental subjects were two-year-old rams and were fed from weaning until slaughter with corn and green hay, with free access to water. Animal were de ned as being unrelated if they had not shared a common ancestor for a minimum of three generations. Within 30 minutes of slaughter, samples of tail adipose tissue were collected from all sheep, snap-frozen, and stored at −80 °C until analysis.
RNA isolation, library preparation, and sequencing A mortar and pestle were used to crush frozen tail fat tissue samples, after which TRIzol (Invitrogen, USA) were used toextract RNA from these samples based on provided directions. A Bioanalyzer 2100 instrument (Agilent Technologies, CA, USA) with an RNA 6000 Nano kit (Agilent, UK) and agarose gel electrophoresis were then used to analyze RNA quantity and quality with high-quality RNA subsequently being used for library preparation. Brie y, libraries preparation for mRNA sequencing was conducted with a standard Illumina standard kit as per the TruSeq RNA SamplePrep Guide (Illumina). Sample mRNA was isolated using oligo-dT beads, after which thermal fragmentation was conducted. Next, random primers and a cDNA Synthesis Kit (Invitrogen) were used to prepare cDNA, which was then converted to yield double-stranded cDNA fragments. All library preparation and sequencing were conducted as per the protocols of the Illumina NextSeq500 instrument (Illumina, USA) in order to generate 75 bp pair-end reads. Small RNA libraries for each of these experimental animals were conducted as above, instead utilizing an Illumina Truseq Small RNA Preparation kit according to the provided instructions. Personalbio (Shanghai, China) performed all sequencing analyses.
The Illumina's Sequencing Control Studio software version (SCS v2.8) were used to obtain raw sequencing data, and base-calling was conducted via the Illumina Real-Time Analysis version 1.8.70 (RTA v1.8.70) software. Extracted sequencing reads were then analyzed.
Mapping, assembly, and annotation Generated mRNA reads were rst quality-checked and converted into FastQC sequence les, after which Trim Galore was utilized with default settings to remove all sequencing adaptors, poly-A tails, and poly-T tails such that only paired-end reads with both pairs > 100 bp were retained for subsequent analysis. Cleaned reads were then mapped against the ovine genome (Oar_v3.1) and the annotation database ENSEMBL with the TopHat v2.0.1 with Bowtie2 software [17,18]. Cu inks v2.0.2 was used to assemble and quantify transcripts with a minimum alignment counts per locus [18,19].
Cleaned unique small RNA reads were obtained by rst conducting quality control with the fastx toolkit (v0.0.13.2), after which Mirdeep2 was used map these reads to the ovine genome assembly (Oar_v3.1, released Sept. 20, 2012) [20]. GenBank and Rfam were then used to exclude corresponding non-coding RNA sequences, and remaining small RNA sequences were searched against the mature miRNAs in miRBase v. 17.0 [21] to identify known miRNAs. Only perfectly mached miRNAs that identical to known precursor or mature RNA sequences within miRBase were annotated as candidate miRNAs.
Differentially expressed gene identi cation Gene expression was assessed in terms of reads per kilobase per million mapped reads (RPKM), and was adjusted via a scaling normalization approach [22], with only those genes that had an RPKM > 1 or < -1 in a minimum of one sequenced sample being considered for downstream analyses. DEG identi cation was conducted with DEseq [23] and Cuffdiff [24], with P values being adjusted to control a false discovery rate (FDR) at 0.05 [25]. DEGs were identi ed as gene with a minimum fold-change of 2 and an adjusted P < 0.05, and were subjected to hierarchical clustering using the heatmap function in the R Bioconductor package [26].
Identi cation of differentially expressed miRNAs and their target genes DEMs were identi ed using the R DESeq package [23], with normalized miRNA expression level being used to compare differential miRNAs expression levels between the generated libraries. Those miRNAs exhibiting a normalized read count of < 5 transcripts per million (TPM) in all libraries were not included in downstream analysis. DEMs were identi ed as those miRNAs with a greater than two-fold change in expression between groups that were signi cant as per the methods de ned by Audic et al. (1997) method [27]. Putative DEM target genes were then identi ed using miRanda [28].
GO and KEGG analysis DEG and DEM target gene enrichment in speci c gene ontology (GO) categories was assessed using WEGO (Web GO) to assess enrichment for speci c cellular component, molecular functions, and biological processes [29,30], with numbers of genes associated with speci c GO terms [31] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [32,33] being quanti ed. The signi cance threshold for GO and KEGG enrichment analyses was P < 0.05.

Integrative DGG and DEM target genes analysis
Correlation between the expression of DEGs and DEMs were conducted in order to analyze four possible directional relationships (both upregulated, both downregulated, upregulated DEMs and downregulated DEGs, and downregulated DEMs and upregulated DEGs). For individual DEMs, analyses were conducted to determine whether numbers of differentially regulated target target genes were higher than predicted by chance alone using lists of both upregulated and downregulated DEGs. DEM and DEG datasets were then integrated, with only DEGs that were negatively correlated with DEMs being considered as targets for downstream analyses.

RNA-sequencing analysis results
After successfully constructing six transcriptome libraries for mRNA expression pro ling. We obtained 36.61 and 33.86 million raw sequence reads for each LTH and STH sample, respectively (  (Table 1).
An average of 16 million raw reads per library was obtained following miRNA library sequencing, yielding 12,248,962, 15,829,550, 16,172,015, 15,480,083, 12,823,610, and 12,314,699 clean reads that were successfully mapped to the sheep genome for the L1, L2, L3, S1, S2, and S3 libraries, respectively ( Table 2). Of these cleaned reads, 27,123 and 27,072 from LTH and STH sheep, respectively, aligned to unique miRNA sequences ( Table 2). Total and unique RNA frequency distributions are shown in Figure 1, revealing highly diverse size distribution mapping to the sheep genome (Fig. 2). Most reads were between 18 to 36 nucleotides-long reads, with 22 nucleotides-long reads being the most common, followed by 21, 23, and 32 nucleotides-long (Fig. 3), consistent with the typical length of miRNAs. Proportions of small non-coding RNAs expression levels in these samples are shown in Table 3. Of these RNAs, miRNAs were the most abundantly expressed in libraries from both LTH and STH sheep (mean = 72.8%; Fig. 1), whereas other RNAs including ribosomal RNAs (rRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs), and transfer RNAs (tRNAs) comprised 27.2% of the total expressed reads. Repeatassociated RNAs made up just 3% of these unique tags, while 7% of these tags did not map to any small non-coding RNA databases or to the sheep genome.
While miRNA do not directly encode proteins, their expression levels can offer functional insights into their role as post-transcriptional regulators. As miRNA functionality is also dependent upon spatiotemporal expression of target mRNAs, we therefore conducted integrated analyses of mRNA and miRNA pro les for LTH and STH sheep. In total, we identi ed 132 miRNAs in the tail fat tissue of these sheep (Table 4), with an average of 130 known miRNAs per sample (range: 125-132). Maxium miRNAs numbers were detected in sheep S2, while the lowest number of miRNAs was detected in the sample from sheep L2.

Analysis of differentially expressed genes between STH and LTH sheep
We next utilized the Cuffdiff software to identify DEGs when comparing STH and LTH sheep samples [19]. In total, we detected 20,137 mRNAs that were expressed in these samples (RPKM>1), including the majority of annotated sheep reference genes (Additional le 1: Table S1). We then assessed correlation between gene expression pro les in different samples, revealing generally strong Pearson correlation coe cients (ranging from 0.724 -0.916; Fig. 4). When selecting DEGs, STH sheep were treated as the control group and LTH sheep as the experimental group. Clustering analyses (Fig. 5A), volcano plots ( Fig.   6) were used to identify DEGs that were up-and down-regulated in LTH sheep relative to STH sheep. In total, we detected 521 signi cant DEGs in tail adipose tissue samples from these sheep (Additional le 2: Table S2), of which 237 were up-regulated and 284 were down-regulated. Several to DEGs associated with fatty acid synthesis, including ADIRF, HSD17B12, LPL, APOBR, INSIG1, THRSP, ACSL5, FAAH, ACSS2, APOA1, ACLY, and ACSM3, were upregulated in LTH sheep, potentially explaining this phenotype.
Analysis of differentially expressed miRNAs and their target genes Next, we quanti ed differential miRNA expression patterns in the tails of LTH and STH sheep in terms of reads per million (RPM) values. All miRNAs with an RPM>10 were annotated based upon mature miRNA sequences recorded in miRBase (release 18.0), leading to the identi cation of 143 mature miRNAs in these six tissue samples (Additional le 3: Table S3). All miRNAs with < 2 nucleotide substitutions outside of the seed region were considered to represent a single miRNA family for the purposes of calculating differential expression. When we assessed tail fat miRNA expression pro les for LTH and STH sheep (Fig.  5B), we detected 132 known miRNAs. Those miRNAs with an average PRM>10 and a fold change between groups of >1 or < -1 were then identi ed as DEMs. Relative to LTH samples, STH samples exhibited 14 DEMs (6 upregulated and 8 downregulated) (Additional le 4: Table S4). The putative targets of these DEMs were then identi ed using miRanda program [28], leading to the identi cation of 2,409 possible target genes (Additional le 5: Table S5). As such, these identi ed DEMs were associated multiple targets and a range of regulatory modules, with an average of 172 target genes per DEM (range: 14-357). Identi ed DEM target genes include those associated with carbohydrate metabolism, the formation of anatomical structures, morphogenesis, lipid metabolism, kinase activity, and immune and in ammatory responses. Identi ed DEM target genes related to lipid metabolism are ACSL4, FTO, FGF8, IGF2, GNPDA2, LIPG, PRKAA2, ELOVL7, SOAT2, and SIRT1.

GOanalysis and functional annotation of DEGs and DEM target genes
A GO analysis of DEGs (RPKM<0.01or |log2Ratio| ≥ 1) was next conducted, with the resultant enriched biological process, cellular component, and molecular function GO terms being summarized in Table S6 (Additional le 6). With respect to biological processes, these genes were highly enriched for genes associated with translation (P=3.29E-05) and lipid metabolic processes (P=9.97e-05). These DEGs were enriched for cellular component GO terms linked to the cytoplasm, extracellular region, extracellular space, and organelles, while enriched molecular function included various structural molecules or functional activities. Overall, the most signi cantly enriched GO terms included extracellular region (P=8.60E-11), external encapsulating structure (P=4.75E-08), and extracellular space (P= 4.1E-06).
To better understand the functional implications of DEM-mediated the target gene regulation, we focused on DEM target genes that were negatively correlated with corresponding DEMs and subjected these genes to a GO analysis (Additional le 7; Table S7). These target genes were signi cantly enriched in cell death processes (P=4.3E-09), and were also associated with carbohydrate metabolism, ligase activity, lipid metabolic processes, protein maturation, and response to stress.

KEGG pathway analysis of DEGs and DEM target genes
Next, we conducted a KEGG pathway analysis of DEGs and DEM target genes identi ed through the above pairwise comparisons, revealing them to cluster primarily into the metabolism, cellular processes, environmental information processing, human diseases, organismal system, and genetic information processing pathways (Additional le 8: Table S8, Additional le 9: Table S9). Many of the identi ed DEGs were particularly enriched in the metabolism and human diseases pathways, when comparing LTH sheep samples with those from STH sheep, the top signi cantly enriched pathways included translation (P=3.4E-07), lipid metabolism (P=4.91E-04), infectious diseases (P=7.1E-09), cancers (P=6.0E-06), carbohydrate metabolism (P=8.4E-06), cardiovascular diseases (P=4.2E-07), and xenobiotic biodegradation and metabolism (P=1.6E-07) (Additional le 8: Table S8). Subsequent KEGG enrichment analysis of the 2,409 identi ed DEM target genes revealed them to be enriched in three metabolismrelated pathways (P < 0.05), including biosynthesis of other secondary metabolites, folding, sorting and degradation, cardiovascular diseases, carbohydrate metabolism, and immune diseases (Additional le 9: Table S9).
In total, these 14 DEMs were associated with 2,409 target genes. We then conducted an integrated analysis of these DEMs, the expression of their target genes, DEGs, and correlations between these expression pro les in sheep tail adipose tissues. We identi ed 65 total miRNA-mRNA interaction pairs through these analyses DEMs and DEGs (Additional le 10: Table S10).

Discussion
Genetics play a major role in determining patterns of fat distribution in animals. LTH sheep are an indigenous sheep breed in China that exhibit a desirable fatty large-tailed phenotype. Such tail fat deposition may be linked to differences in metabolic activities, adipocyte size, or developmental factors. Indeed, tail adipose tissue develops at later those in renal or cardiac adipose reserves [34]. Adipose tissue expansion is closely linked to the e cient storage and release of fatty acids within adipocytes, but the genetic genetic basis for the fatty tail phenotype observed in LTH sheep remains to be de ned. We therefor conducted an RNA sequencing-based transcriptomic analysis of mRNA and miRNA expression patterns in adipose tissues from the tail of six LTH and STH sheep.
First, we successfully sequenced six transcriptomeic libraries and obtained 36.61 and 33.86 million raw sequence reads for our different sample types (Table 1), with over 94% of these raw reads being uniquely aligned to the Oar_v3.1 sheep reference genome, and with very high mapping rates. Based on these metrics, our results were more reliable than those from prior studies of Kazak and Tibetan sheep [14] or Han and Dorset sheep [16]. We therefor concluded that our data provided su cient coverage to permit downstream analyses, as an average of 30 million reads has been shown to provide >90% coverage of annotated genes in yeast [35], chicken [36], and cattle [37]. These result were in line with those of prior studies of other sheep breeds [14,16]. In miRNAs, an average of 16 million raw reads per library was obtained ( Table 2). Raw reads after quality control, and these cleaned reads aligned to unique miRNA sequences (Table 2). RNA size distributions and mapping to the sheep genome are shown in Figure 1 and 2. Most reads were between 18 to 36 nucleotides-long reads (Fig. 3). Proportions of small non-coding RNAs expression levels in these samples are shown in Table 3. Of these RNAs, miRNAs were the most abundantly expressed in libraries from both LTH and STH sheep. While miRNA do not directly encode proteins, their expression levels can offer functional insights into their role as post-transcriptional regulators. We therefore identi ed 132 miRNAs in the tail fat tissue of these sheep (Table 4).  previously identi ed 3,132 miRNAs in the adipose tissue of STH and Dorset sheep via RNAsequencing [15]. In contrast, we identi ed signi cantly fewer miRNA in the present study (Additional le 3: Table S3), potentially due to differences in breeds, sequencing depth, and/or ltering criteria.
When selecting DEGs, STH sheep were treated as the control group and LTH sheep as the experimental group. In total, we detected 20,137 mRNAs that were expressed in these samples (Additional le 1: Table   S1), and of these mRNAs, 521 signi cant DEGs were found in tail adipose tissue samples from these sheep (Additional le 2: Table S2), of which 237 were up-regulated and 284 were down-regulated. Clustering analyses (Fig. 5A) and volcano plots (Fig. 6) were used to identify DEGs that were up-and down-regulated in LTH sheep relative to STH sheep. In a prior study of Kazak and Tibetan sheep, authors identi ed 646 DEGs including 280 and 366 that were up-and down-regulated, respectively [14]. Miao et al.
(2015) similarly identi ed 602 DEGs when comparing samples from STH and Dorset sheep [16]. As such, we detected fewer DEGs than in these prior studies (Additional le 2: Table S2), suggesting that distinct regulatory mechanisms govern adipogenesis in these different breeds.
In our identi ed DEGs, such as ADIRF (Adipogenesis Regulatory Factor) encodes a protein that is closely linked to adipocytic differentiation and development, controlling transcription activity early during the preadipocyte differentiation process [38,39]. HSD17B12 (Hydroxysteroid 17-Beta dehydrogenase 12) has been associated with suppressed adipocyte lipid accumulation, consistent with its predicted effect predicted in model system [40]. One key step in the energy metabolism process is the hydrolysis of triacylglycerol-rich lipoproteins by lipoprotein lipase in the vascular endothelium, resulting in the release of fatty acids that can subsequently be metabolized or stored for future use. APOBR (Apolipoprotein B receptor) encodes a receptor expressed on macrophage that controls cellular fat and vitamin uptake [41,42]. In normal weight humans, a single high fat meal increased both APOBR expression and lipid uptake in monocytes [43]. High blood lipid levels differentially regulate APOBR expression in human postprandial monocytes and macrophages and lead to foam cell formation [44]. INSIG1(Insulin-induced gene 1) is a transmembrane ER-associated protein that retains sterol regulatory element-binding transcription factors (SREBFs) in the ER, thus preventing them from undergoing proteolytic activation within the Golgi apparatus and thereby controlling lipid metabolism [45,46], INSIG1 also in uences human obesity-related hypertriglyceridemia [47].THRSP (Thyroid hormone-inducible hepatic protein) plays a role in regulating fatty acid synthesis-related gene expression [48], and its bovine homolog has been linked to muscle fatty acid composition in cattle [49]. In mice, ACSL5 (Acyl-CoA synthetases 5) regulated systemic energy metabolism such that the knockout of this gene can prevent the onset of obesity and insulin resistance [50]. DGAT2 (Diacylgycerol acyltransferase 2) is an essential catalyst of triglyceride biosynthesis, which is important given that excessive triglyceride accumulation within adipose tissues is a key facet of obesity [51]. Overexpression of FAAH (Fatty acid amide hydrolase) can suppress DGAT2 expression and triglyceride synthesis, and these genes may interact with one another to in uence adiposity [51]. ACSS2 (Acetyl-CoA synthetase short-chain family member 2) codes for a cytosolic enzyme that catalyzes acetate activation for use in the context of lipid synthesis and energy production. Vysochan et al. (2017) found that ACSS2 was able to convert glucose-derived carbon into acetate for the generation of cytosolic acetyl-CoA during lipid synthesis, which they found to be a key process in the replication of HCMV and in virus-induced lipogenesis [52]. ApoA-I/HDL levels have also been linked to decreased body weight and enhanced insulin sensitivity [53,54], such that ApoA-I-de cient mice exhibited increased body weight and body fat levels as well as decreased weight loss during caloric restriction, at least in part due to impaired lipolysis [55,56]. Berton et al. (2016) determined that ACSM3 (Acyl-CoA synthetase medium-chain family member 3) was liked to lipid metabolism and fatty acid composition in Nellore cattle, thereby impacting key meat quality traits [57]. Microarray analyses conducted in mice have additionally revealed that genes associated with lipolysis, fatty acid metabolism, mitochondrial energy transduction, oxidation-reduction, insulin sensitivity, and skeletal system development are downregulated in animals fed a high fat diet, whereas genes related to extracellular matrix remodeling and in ammation, including ACSM, are upregulated in these animals [58]. ACL (ATP citrate lyase) function in the cytosol wherein it converts mitochondria-derived citrate into oxaloacetate and acetyl-CoA, which can then be leveraged to facilitate lipid synthesis and acetylation [59]. Zaidi et al. (2012) also determined that a lack of ACLY bolstered the ACSS2-dependent synthesis of lipids [60]. Adipocytes lacking ACLY expression in vivo accumulate lipids, produce higher levels of acetate-derived acetyl-CoA and malonyl-CoA, and exhibit distinct patterns of fatty acid content and synthesis [61]. These genes were associated with lipid metabolism and fatty acid synthesis, potentially explaining this phenotype. This suggests that DEGs linked to lipid synthesis and triglyceride hydrolysis may mediate the enanced accumulation of lipids within the tails of LTH sheep, highlighting these genes as important targets for further study.
We quanti ed differential miRNA expression patterns in the tails of LTH and STH sheep in terms of reads per million (RPM) values. 143 mature miRNAs were identi ed in these six tissue samples (Additional le 3: Table S3). We assessed tail fat miRNA expression pro les for LTH and STH sheep (Fig.  5B), and identi ed 132 known miRNAs. Relative to LTH samples, STH samples exhibited 14 DEMs (6 upregulated and 8 downregulated) (Additional le 4: Table S4), with 2,409 putative targets having been identi ed (Additional le 5: Table S5). We found that identi ed DEMs played a range of roles in controlling tail fat deposition, development, and metabolism. For example, the Let-7 family of miRNAs have been shown to control glucose metabolism and blood glucose levels in murine loss-of-function experiments, at least in part owing to its ability to regulate muscle insulin signaling [62]. We found that oar-let-7i was downregulated in the present study, suggesting that it may control lipid metabolism and fat deposition in sheep tails. Additionally, miR-133 antagonism in the context of muscle regeneration has been shown to enhance respiratory uncoupling, thermogenesis, and glucose uptake within muscle, in addition to improving overall glucose tolerance, energy expenditure, and diet-induced obesity resistance [63]. Levels of mir-133 are decreased in mice that are exposed to cold, leading to de novo satellite cellderived brown adipocyte generation [63]. We determined that oar-mir-133 was upregulated, highlighting it as another putative regulator of ovine tail fat accumulation. Stable miR-376a and miR-376c overexpression have been shown to slow growth and migratory activity in vitro, and both of these miRNAs are predicted to suppress IGF1R. Such suppression has been con rmed via luciferase reporter assay in melanoma cells [64]. miR-411 is a miR-379 family member encoded in the DLK-DIO3 region of human chromosome 14 [65]. Harafuji et al. (2013) determined that miR-411 exhibits differential expressed in FSHD myoblasts, suggesting it may function as a controller of myogenesis [66]. miR-221 and miR-222 are also associated with myogenesis and are necessary for fully myocyte differentiation [67]. We found oar-miR-221 to be differentially regulated in the present study, suggesting that it may serve as a regulator of ovine fat deposition. Lee et al. (2014) found that miR-543 and miR-590-3p directly suppressed AIMP3/p18 translation, thereby impairing the onset of senescence [68]. Formosa et al. (2014) further found that miR-154, miR-376a, miR-376c, miR-381, and miR-495 controlled metastatic prostate regulate cancer cells proliferation, apoptosis, migration, and invasion [69]. Additional research has identi ed miR-495 as an inhibitor of bone regeneration, gastric cancer cell migration, and invasion, potentially owing to its ability to target high-mobility group AT-Hook 2 [70,71]. We identi ed oar-miR-150, oar-miR-30a-3p, oar-miR-543-3p, oar-miR-412-3p, and oar-miR-154a-3p as DEMs in the present study, suggesting that they may be important regulators of metabolic activity and tail fat deposition in sheep. Identi ed DEM target genes included genes associated with carbohydrate metabolism, the formation of anatomical structures, morphogenesis, lipid metabolism, kinase activity, and immune/in ammatory responses. Overall, these results suggested that the identi ed DEMs may play central roles in the regulation of metabolism-associated gene expression in the context of sheep tail fat deposition, with both DEGs and DEMs being closely linked to metabolic activities including lipid metabolism.
We additionally leveraged GO enrichment analyses to identify certain genes that were important DEM regulatory targets in the context of lipid metabolism. Of the identi ed DEMs, miR-133 was predicted to target the greatest number of lipid metabolism-related genes, including ACSL4, FTO, FGF8, IGF2, and GNPDA2. ACSL4 is an essential regulator of fatty acid metabolism, and the expression of ACSL4 is controlled by PPARδ-mediated regulatory [72]. Smemo et al. (2014) established IRX3 as a functional longrange target of obesity-associated FTO variants and a key regulator of overall body composition [73]. Lopez-Sanchez et al. (2015) also determined that miR-133 is an essential mediator of crosstalk between the Fgf8 and Bmp2 signaling pathway owing to its ability to regulate Fgf8/ERK signaling in the context of cardiac induction [74], in line with our identi ition of Fgf8 as a miR-133 target genes. IGF2 is a DEM target gene that has likely been subjected to recent selective pressure, given that most domestic European pigs selected for lean growth carry a particular IGF2 allele that induces decreased fat deposition and increased muscle growth [75]. The GNPDA2 (glucosamine-6-phosphate deaminase 2) gene is an allosteric enzyme responsible for catalyzing reversible D-glucosamine-6-phosphate conversion into D-fructose-6-phosphate and ammonium. GNPDA2 gene variants have previously been linked to obesity susceptibility and body mass index [76], and human GWAS studies have supported these links [77]. LIPG (Endothelial lipase) is a phospholipase that can decrease plasma HDL cholesterol levels [78]. ELOVL7 (Elongation of very-long chain fatty acids protein 7) controls the elongation of C18:0, C20:0, and other very-long-chain fatty acids [79], with in vitro analyses having shown this gene to speci cally control saturated very-long-chain fatty acids elongation [80]. Prior studies of pigs have identi ed ELOVL7 as a candidate gene related fatty acid composition in porcine muscle and abdominal fat tissues [81], suggesting that it also warrants research to evaluate its impact on long-chain fatty acids in sheep. PRKAA2 (Protein Kinase AMP-Activated Catalytic Subunit Alpha 2) encodes a catalytic subunit of the AMP-activated protein kinase (AMPK), which is a central regulator of cellular energy status, and involved in the glucose and lipid metabolism [82]. Lee et al. (2007) have investigated the porcine PRKAA2 gene as a positional candidate regulator of intramuscular fat and backfat thickness traits in pig chromosome 6 [83]. SOAT2 (Sterol O-acyltransferase 2), also known as ACAT2, encodes an enzyme that is primarily expressed in the intestine and liver and that controls cholesterol esteri cation, beta-oxidation, and lipid metabolism [84]. It serves as a key regulator of lipid and cholesterol metabolism in humans and cattle, and signi cantly ACAT2 variability has been detected in samples of human liver tissue, indicating that it is associated with over half of ACAT activity in most humans [84]. In one study of broiler chickens, hepatic ACAT2 expression was shown to regulate abdominal fat accumulation of [85]. Correlations between cholesterol and ACAT2 have been reported in sh [86], and in human have been linked to increased lean mass [87]. ACAT2 has also been shown to regulate enzymes associated with lipogenesis and lipolysis in beef cattle and pigs [88,89]. SIRT1 (sirtuin 1) regulates lipid and glucose metabolism to control overall energy metabolism, and genetic research has been conducted to assess the impact of SIRT1 SNPs on adiposity. In one casecontrol association study conducted in Belgium of 1,068 obese patients and 313 lean controls, SIRT1 variants were linked to an increased risk of obesity and were correlated with visceral obesity-related parameters in obese males [90]. Both SNPs and overall SIRT1 expression levels have also been linked to severe obesity [91]. Together, our results suggest that these DEMs and their target genes are likely to play ket roles in regulating sheep tail fat deposition and lipid metabolism.
Many of the genes enriched in the identi ed pathway were associated with lipid synthesis and metabolisms, with the majority also being associated with the fatty acid metabolism pathway. SIRT1 has previously been identi ed as a key regulator of cellular differentiation, energy metabolism, and resistance to stress [92]. In total, we additionally identi ed 28 genes that were linked to adipogenesis and lipid metabolism in prior studies. Many of these genes were downregulated in tail fat tissue samples from LTH sheep relative to those from STH sheep, consistent with the hypothesis that lipid metabolic activities are decreased in tail adipocytes of LTH sheep. In addition, we found that certain tail fat deposition-related genes (FASN, FABP4, LPL, THRSP, and DGAT1) were also downregulated in the tail fat of these animals. These patterns of gene expression were distinct from those in similar prior studies, and emphasize the importance of conducting further research to better understand the degree to which these differences are attributable to breed-speci c differences. Additionally, it may be bene cial to evaluate the functional relevance of tail fat deposition-related gene expression in different sheep breed by grouping tail fat samples according to corresponding tail weight indices.
As detailed above, when we conducted an integrated analysis of these DEMs, the expression of their target genes, DEGs, and correlations between these expression pro les in sheep tail adipose tissues. We identi ed 65 total miRNA-mRNA interaction pairs through these analyses DEMs and DEGs (Additional le 10: Table S10). Given the large number of target genes identi ed, it is likely that this small set of DEMs can regulate the expression of many important genes in the context of tail fat deposition. miRNAs are short non-coding RNAs that post-transcriptionally regulate gene expression. They generally bind to complementary sequence in the 3'-untranslated region (UTR) of their target mRNAs and repress protein production via translational silencing and mRNA destabilization [93]. In addition, long noncoding RNAs (lncRNAs) can additionally function as transcriptional regulators by binging to histone-modifying complexes, DNA binding proteins, and even to RNA polymerase II, thereby modulatin gene expression [94]. While we identi ed 14 relevant DEMs and 521 DEGs in this study, relatively few interaction pairs were identi ed based on these sequencing data. Differences obtained when using different DEG lists suggest that the miRNA control mRNA expression in a carefully regulated manner.

Conclusions
In summary, in the present study, we conducted an integrated analysis of mRNA and miRNA expression patterns in tail adipose tissue from LTH and STH sheep, leading to the identi cation of DEGs and DEMs that were highly expressed in LTH tissue samples relative to STH samples. These DEGs and DEMs were predicted to play key functional and regulatory roles in the regulation of deposition and metabolic activation in the tails of LTH sheep. Overall, these data provide a foundation for further studies of the biological roles of these DEGs, DEMs, and DEM target genes in the context of ovine tail fat deposition, highlighting potential mechanisms governing this fatty large-tailed phenotypic trait.

Declarations
The experiment was approved by the Ministry of Science and Technology of China and the Ethics Committee of Shangqiu Normal University. All animal procedures were conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of China (Guidelines on Ethical Treatment of Experimental Animals (2006) No. 398).

Availability of data and materials
The datasets used and analysed during this study are available on request from the corresponding author. Figure 1 Frequency of small RNAs in tail fat tissues of the LTH and STH sheep breeds. The proportion of expressed adipose tissues small RNAs identical to LTH and STH sheep breed: i) the known miRNAs are shaded light green, ii) the novel miRNAs are shaded red, iii) repeat-associated RNA sequences are shaded blue (including rRNA, snoRNAs, and snRNAs), iv) tRNAs are shaded yellow, and v) the proportion of mapped but unclassi ed sequences is shaded brown.