Abstract

   microRNAs are a class of important non-coding RNAs, which can
   participate in the regulation of biological processes. In recent years,
   miRNA has been widely studied not only in humans and mice, but also in
   animal husbandry. However, compared with other livestock and poultry
   breeds, the study of miRNA in subcutaneous adipose tissue of sheep is
   not comprehensive. Transcriptome analysis of miRNAs in subcutaneous
   adipose tissue of Duolang sheep, and Small Tail Han sheep was performed
   using RNA-Seq technology. Differentially expressed miRNAs were screened
   between different breeds. Target genes were predicted, and then the
   joint analysis of candidate genes were conducted based on Gene Ontology
   (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment.
   Finally, the RNA-Seq data were verified by real-time quantitative
   polymerase chain reaction (qRT-PCR). Herein, we identified 38
   differentially expressed miRNAs (9 novel miRNAs and 29 known miRNAs).
   In addition, a total of 854 target genes were predicted by miRanda
   software. GO and KEGG pathway analysis demonstrated that regulation of
   lipolysis in adipocytes plays a key role in the deposition of
   subcutaneous adipose tissue in Duolang sheep and Small Tail Han sheep.
   The miRNAs might regulate fat deposits by regulating genes involved in
   regulation of lipolysis in adipocytes. Specifically, NC_ 040278.1_
   37602, oar-mir-493-3p, NC_ 040278.1_ 37521 and NC_ 040255.1_ 11627
   might target PTGS2, AKT2, AKT3, and PIK3CA, respectively, and then play
   critical regulatory role. In conclusion, all the results provide a good
   idea for further revealing the mechanism of subcutaneous adipose tissue
   deposition and improving the meat production performance of sheep, and
   lay a foundation for promoting the development of animal husbandry.

   Keywords: sheep, subcutaneous fat, microRNAs, high-throughput
   sequencing, fat deposition

Introduction

   miRNA is a kind of short-chain endogenous non-coding RNA with a length
   of 18–25 nt, which is highly conservative and widely distributed in
   animals, plants and viruses ([31]1, [32]2). With the discovery of lin-4
   and let-7 in Caenorhabditis elegans ([33]3, [34]4), the research of
   miRNA has gradually entered people's vision. Although miRNAs cannot
   encode proteins, they still control various biological processes and
   are key regulators of development and cellular homeostasis ([35]5). In
   many biological processes, miRNAs regulate gene expression by
   recognizing homologous sequences and interfering with transcriptional,
   translational or epigenetic processes ([36]6). In the past few years,
   the research of miRNA in cancer has exploded ([37]7), including
   diseases such as hepatocellular carcinoma ([38]8), breast cancer
   ([39]9), and gastric cancer ([40]10). Not only that, miRNA in type 2
   diabetes ([41]11), hypertension ([42]12) and atherosclerosis ([43]13)
   is also very common. Obesity has become a prevalent health problem as
   it is closely associated with many diseases.

   Obesity is due to the proliferation, differentiation or enlargement of
   adipocytes, and miRNAs have both promoting and inhibiting effects on
   adipocyte differentiation ([44]14). LPL is a direct target gene of
   miR-152 during preadipocyte differentiation. MiR-152 can inhibit the
   proliferation of 3T3-L1 preadipocytes and promote the differentiation
   of 3T3-L1 preadipocytes by negatively regulating LPL ([45]15). The
   study found that miR-244 regulates the adipogenic differentiation of
   bovine preadipocytes by targeting LPL. When miR-224 was overexpressed,
   the mRNAs expression levels of the adipogenesis-related markers PPARγ,
   FASN, C/EBPα, C/EBPβ and PLIN1 decreased, whereas the opposite effect
   was produced when miR-224 was inhibited ([46]16). The regulation of
   lipid metabolism by miRNA is a complex network regulation system, which
   mainly affects lipid metabolism by regulating genes related to lipid
   synthesis, oxidation, and transport ([47]17, [48]18). For example,
   miR-33, as one of the widely studied miRNAs, is involved in the
   regulation of multiple lipid metabolism links. The expression of genes
   related to fatty acid β-oxidation can be inhibited, and the
   biosynthesis of high-density lipoprotein and the outflow of cholesterol
   can also be inhibited ([49]19, [50]20). As the first miRNA found to be
   involved in the regulation of lipid metabolism, miR-122 has been deeply
   studied. MiR-122 is most abundant in the liver and plays an important
   role in mainly regulating lipid synthesis and oxidation processes,
   especially the accumulation of triglycerides. MiR-122 plays an
   important regulatory role in the translation of YY1 and FXR. Mir-122
   upregulates FXR expression by targeting the 3'UTR of YY1 mRNA upstream,
   suppressing triglyceride levels in hepatocytes, and FXR is also
   involved in the regulation of lipid metabolism disorders and insulin
   resistance ([51]21). Another study showed that the expression of
   miR-122 was inhibited by binding to the 3'UTR of Sirt1 in non-alcoholic
   fatty liver disease, and down-regulation of miR-122 inhibited
   adipogenesis genes, but activated the AMPK signaling pathway, further
   inhibiting hepatic lipogenesis and triglyceride secretion ([52]22).
   Studies have found that miR-122 can also promote cholesterol synthesis
   ([53]23, [54]24). The above studies show that miRNA plays a key
   regulatory role in the process of fat deposition.

   As one of the main meat livestock and poultry resources in the world,
   sheep has warm meat, high protein content and low fat and cholesterol
   levels. Compared with pork, the meat quality is more delicate, and the
   content of essential amino acids is also higher than that of pigs,
   chickens and cattle ([55]25, [56]26). Adipose tissue is an important
   factor affecting meat quality, mainly including the effects on sensory,
   flavor and tenderness ([57]27, [58]28). Moreover, the back-fat of sheep
   is an important indicators of meat yield and hot carcass composition
   ([59]29). And subcutaneous adipose tissue accounts for the highest
   proportion of total fat content in animals ([60]30), which has the
   function of protecting animals and storing energy. Duolang sheep and
   Small Tail Han sheep belong to high-quality local breeds in China, and
   there are differences in fat deposition between the two breeds. Duolang
   sheep have large fat buttocks, strong meat production capacity, fresh
   and juicy meat, and belongs to both meat and fat sheep, and Duolang
   sheep also have the characteristics of rapid growth and development.
   Small Tail Han sheep are short and thin-tailed sheep, with strong
   environmental adaptability and reproductive ability, fast growth rate
   and stable genetic performance, but the meat body size is not obvious
   and the carcass meat production rate is low. At present, living
   standards have been greatly improved, and people's requirements for
   meat quality and quantity have also increased. Therefore, studying the
   molecular mechanism of sheep fat deposition can improve the quality of
   meat products to meet consumer demand. At the same time, adipose tissue
   is an important organ of energy metabolism, and excessive fat
   deposition can lead to obesity and the occurrence of a series of
   metabolic syndromes ([61]31). Therefore, studying the mechanism of fat
   deposition can not only improve meat quality and promote the
   development of animal husbandry, but also prevent or treat a series of
   diseases affecting human health caused by excessive fat deposition. We
   selected the subcutaneous adipose tissue of Duolang sheep and Small
   Tail Han sheep as experimental materials, analyzed the gene expression
   profile of subcutaneous adipose tissue by transcriptome sequencing and
   bioinformatics methods, screened and identified the key candidate genes
   related to lipid metabolism and adipogenic differentiation and explored
   the molecular mechanism related to fat deposition.

Materials and Methods

Sample Collection and RNA Isolation

   All work in this study was approved by the Animal Welfare and Ethics
   Committee of Beijing Institute of Animal Sciences, Chinese Academy of
   Agricultural Sciences (No. IAS2019-82). In order to detect the
   expression profile of miRNAs in sheep subcutaneous fat, we collected
   subcutaneous adipose tissue samples from 3 adult female Duolang sheep
   (D-PF-1, 2, 3) and 3 adult female Small Tail Han sheep (X-PF-1, 2, 3).
   All sheep were raised in the same conditions with free to drink and eat
   under natural light, and their diet meets the current nutritional
   needs. They were in good physical condition, aged 2 years old and the
   mean bodyweight of sheep were 50 ± 3 kg. We collected subcutaneous
   adipose tissue samples located at the back fat according to the
   agricultural industry standard of the people's Republic of China (NY/T
   3469-2019). In order to reduce the pain of animals, we first stunned
   them with electricity and then slaughtered them. The whole sampling
   process was controlled within 30 min. The collected samples were
   immediately frozen in liquid nitrogen and stored in an environment of
   −80°C before the experiment.

   Total RNA was isolated from each adipose tissue using TRIzol reagent
   (Invitrogen, Invitrogen Life Technologies, Carlsbad, USA), and genomic
   DNA was removed using rDNase I RNase-free (TaKara). The RNA
   concentration and quality were then measured with a 2100 Bioanalyzer
   (Agilent Technologies, Santa Clara, CA, USA) and the ND-2000 (NanoDrop
   Technologies), RIN > 8 was good quality.

Library Preparation and Sequencing

   The RNA library was constructed using the TruSeqTM Small RNA sample
   prep Kit (Invitrogen) according to the instructions. The rRNA in the
   total RNA was first removed, the 3′ end adapter and the 5′ end adapter
   were, respectively, connected with the kit, and then the random primers
   were reversed to 1st cDNA. To enrich the library we performed 11-12 PCR
   cycles. The library was enriched and then purified (6% Novex TBE PAGE
   gel, 1.0 mm, 10 well) and quantified by TBS380 (Picogreen). Bridge PCR
   was performed on cBot to generate clusters. The library was sequenced
   with the Illumina NovaSeq 6000.

Quality Control and Sequence Alignment

   In order to improve the sequencing quality and get clean reads,
   introduced adapter sequences, low-quality reads, sequences with high N
   rate (N stands for indeterminate bases), and sequences that are too
   short will be removed. The clean reads were aligned with the reference
   genome using Bowtie ([62]32) software. The reference genome
   GCF_016772045.1 was obtained in the NCBI database. StringTie ([63]32)
   software was used to splicing the mapped reads.

miRNA Identification and Structural Analysis

   The reads aligned to the reference genome were aligned with the miRNA
   precursor and mature sequences in the miRBase V.22.1 ([64]33) to obtain
   known miRNAs. The Rfam ([65]34) was used to filter ncRNAs and repeats
   sequences such as ribosomal RNA (rRNA), transfer RNA (tRNA), small
   nuclear RNA (snRNA) and small nucleolar RNA (snoRNA), and at that time,
   the types and numbers of these sequences were counted. The sRNAs that
   cannot be aligned with Rfam and miRBase were aligned to the reference
   genome, and the surrounding sequences were intercepted using miRDeep2
   ([66]35) software for secondary structure prediction to identify new
   miRNAs. In the process of developing from precursor to mature miRNA,
   the Dicer restriction site was specific, which makes the first base of
   miRNA mature sequence strongly biased to U. In addition, some
   non-canonical editing sites exist in miRNAs, thereby altering target
   genes. And the MiRME ([67]36) method was used to detect various
   mutation and editing sites of miRNA. Finally, based on the seed
   sequence, the identified known miRNAs and new miRNAs were subjected to
   miRNA family analysis.

Differentially Expressed miRNAs

   In order to facilitate the subsequent analysis of the differential
   expression among the samples, quantitative analysis was performed on
   the expression levels of the samples, respectively. RSEM ([68]37)
   software was used for quantification, and the expression level of each
   miRNA was calculated according to the transcripts per million readsx
   (TPM) method ([69]38). DESeq2 ([70]39) is a statistical analysis based
   on negative binomial distribution, which can be used in experiments
   with biological replicates. Significant differently expressed miRNAs
   were extracted with p-value < 0.05 and |log[2]FC|≥1.

miRNA Target Gene Prediction

   Because miRNAs cannot encode proteins, they function through
   post-transcriptional regulation of their target genes. In animals,
   miRNA relies on the seed sequence (2–8 nt at the 5′ end) to closely
   bind to the 3′ non-coding region of the target gene, thereby inhibiting
   the translation of the target mRNA. In this study, the miRanda ([71]40)
   software was used to predict its target genes, and the parameters were
   set as: Score ≥ 160 and Energy ≤ −20. The predicted target genes were
   visualized by the application software Cytoscape ([72]41). All of the
   mRNA data was obtained from fat collected from the same animals that
   were used for the miRNA analysis at the same relative time (age) and
   using a similar procedure.

GO and KEGG Enrichment Analysis of Differentially Expressed Genes

   GO ([73]42) can be used for functional enrichment analysis on the
   differentially expressed genes. The software Goatools was used to
   perform GO enrichment analysis on the differentially expressed genes,
   so as to obtain functional annotations of the genes, including three
   parts: biological process, molecular function and cellular components.
   The KEGG ([74]43) pathway enrichment analysis was also performed on the
   genes, and the R software was used for the enrichment analysis of
   metabolic pathways and information processing pathways. Using Fisher's
   exact test and Benjamini and Hoceberg (BH) to correct p-value to obtain
   p-adjust, when p-adjust < 0.05, we considered significant enrichment.

Quantitative Real-Time Polymerase Chain Reaction

   Real-time fluorescent quantitative PCR (qRT-PCR) was used to verify the
   reliability of the sequencing results. Five differentially expressed
   miRNAs and eight differentially expressed mRNAs were randomly selected
   for verification. In total, 0.5 μg RNA was taken to synthesize cDNA
   template through GeneAmp^® PCR System 9700 (Applied Biosystems, USA).
   First, RNA, 4 × g DNA wiper Mix and nuclease free H[2]O were reacted in
   GeneAmp^® PCR System 9700 at 42°C for 2 min. Second, 5 × HiScript II Q
   RT SuperMix IIa were added and reacted at 50°C for 15 min, then for 5 s
   at 85°C. And then the reverse transcription reaction mix were dilute ×
   10 in nuclease free H[2]O and maintained at −20°C. The qRT-PCR analysis
   was conducted using LightCycler^® 480 II Real-time PCR Instrument
   (Roche, Swiss). U6 and ACTB were used as miRNA and mRNA internal
   references, respectively. Three biological replicates were employed for