Abstract

   MicroRNAs (miRNAs) are a class of small non-coding RNAs that have been
   shown to play important post-transcriptional regulatory roles in the
   growth and development of skeletal muscle tissues. However, limited
   research into the effect of miRNAs on muscle development in goats has
   been reported. In this study, Liaoning cashmere (LC) goats and Ziwuling
   black (ZB) goats with significant phenotype difference in meat
   production performance were selected and the difference in Longissimus
   dorsi muscle tissue expression profile of miRNAs between the two goat
   breeds was then compared using small RNA sequencing. A total of 1,623
   miRNAs were identified in Longissimus dorsi muscle tissues of the two
   goat breeds, including 410 known caprine miRNAs, 928 known
   species-conserved miRNAs and 285 novel miRNAs. Of these, 1,142 were
   co-expressed in both breeds, while 230 and 251 miRNAs were only
   expressed in LC and ZB goats, respectively. Compared with ZB goats, 24
   up-regulated miRNAs and 135 miRNAs down-regulated were screened in LC
   goats. A miRNA-mRNA interaction network showed that the differentially
   expressed miRNAs would target important functional genes associated
   with muscle development and intramuscular fat deposition. Kyoto
   Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed
   that the target genes of differentially expressed miRNAs were
   significantly enriched in Ras, Rap 1, FoxO, and Hippo signaling
   pathways. This study suggested that these differentially expressed
   miRNAs may be responsible for the phenotype differences in meat
   production performance between the two goat breeds, thereby providing
   an improved understanding of the roles of miRNAs in muscle tissue of
   goats.

   Keywords: microRNA (miRNA), muscle development, intramuscular fat,
   small RNA sequencing, goat

Introduction

   MicroRNAs (miRNAs) are a class of non-coding small RNA molecules (~22
   nucleotides), which are evolutionarily conserved in eukaryotes ([47]1).
   In recent years, miRNAs are well-recognized as negative regulators of
   gene expression at post-transcriptional level, in that they can either
   inhibit translation or promote degradation of mRNA by complementary
   binding to the 3′-untranslated regions (3′-UTR) of the target genes.
   The miRNAs are therefore involved in a wide variety of cell biological
   processes, including proliferation, differentiation, death, and fate
   specification ([48]1, [49]2).

   In modern animal husbandry, skeletal muscle is considered as the most
   economically important tissue of producing-meat livestock. Many studies
   have confirmed that miRNAs played essential roles in the growth and
   development of skeletal muscle. For example, miR-1 and miR-206 have
   been reported to facilitate differentiation and inhibit proliferation
   of bovine skeletal muscle satellite cells by targeting PAX7 ([50]3).
   The over-expression of miR-486 induced skeletal muscle hypertrophy of
   mice via activation of protein kinase B (Akt) ([51]4). The miR-27b
   regulated caprine myogenic proliferation and differentiation of
   skeletal muscle satellite cells by inhibiting the expression of PAX3
   ([52]5).

   Up to now, research into expression profiles of miRNAs in the skeletal
   muscle tissue of domestic animals have mainly been focused on pigs
   ([53]6–[54]9), cattle ([55]10–[56]13), and sheep ([57]14–[58]17). It
   was found from these studies described above that miRNAs were
   differentially expressed in skeletal muscle at different developmental
   periods, or between different breeds. These further demonstrated
   crucial effect of miRNAs on skeletal muscle development.

   In goats, the studies of miRNA expression profiles in skeletal muscle
   tissue have mainly been focused on different development stages. For
   example, Wang et al. identified 336 differentially expressed miRNAs in
   skeletal muscle of Huanghuai goats between fetal stage and 6-month-old
   stage, of which miR-424-5p, miR-29a, miR-129-3p, miR-181b, and miR-181d
   were involved in multiple important pathways related to muscle
   development ([59]18). Guo et al. and Ling et al. also found some
   important differentially expressed miRNAs in skeletal muscle from
   prenatal stages to neonatal stage in Jianzhou Da′er goats and Anhui
   white goats ([60]19, [61]20). However, little is known about the miRNA
   profiles of muscle tissues in other goat breeds, or between different
   goat breeds.

   Ziwuling black (ZB) goats and Liaoning cashmere (LC) goats are both
   indigenous goat breeds in China. There are significant differences in
   meat production performance and muscle nutrients between the two
   breeds. For example, LC goats had higher carcass weight, muscle mass
   and intramuscular fat content, but lower muscle fiber density as well
   as contents of linoleic (C18: 2n-6), 11C, 14C-eicosadienoic acid (C20:
   2n-6), moisture, and crude ash in meat when compared to ZB goats (P <
   0.05) ([62]21). In this context, elucidating the molecular mechanisms
   regulating these phenotypic differences between LC and ZB goats can
   provide insight for improving meat production performance of goats and
   other livestock. In this study, the expression profiles of miRNAs were
   compared in the Longissimus dorsi muscle between LC and ZB goats using
   small RNA sequencing. The differentially expressed miRNAs were screened
   between the two caprine breeds and the roles of miRNAs were also
   uncovered in skeletal muscle development and intramuscular fat
   deposition in goats.

Materials and Methods

Ethics Statement

   All animal procedures in this study were approved by Animal Experiment
   Ethics Committee of Gansu Agricultural University with an approval
   number of GSAU-ETH-AST-2021-028.

Longissimus dorsi Muscle Sample Collection and RNA Extraction

   Ten healthy, 9-month-old male goats (five LC goats and five ZB goats)
   were selected from Yongfeng Goat Breeding Company in Huan County, Gansu
   Province, China. All goats were raised under the same environmental
   conditions and nutrition levels. After being slaughtered, the
   Longissimus dorsi muscle samples from the area between 12th and 13th
   ribs on the left carcass of each goat were collected and then frozen in
   liquid nitrogen immediately until further use. The meat production
   performance, muscle fiber size and intramuscular fat content from these
   LC and ZB goats have been reported by Wang et al. ([63]21) and were
   also presented in [64]Supplementary Material 1.

   Total RNA from Longissimus dorsi muscle samples was extracted using a
   Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). The concentration
   and purity of the RNA extracted were assessed using a Nanodrop 2000
   (Thermo Scientific, MA, USA). Only samples with an RNA concentration
   >80 ng/uL and a purity of 1.80–2.10 were used for the study. The
   Agilent 2100 Bioanalyzer (Agilent, CA, USA) was used to assess RNA
   Integrity Number (RIN) of samples. Only RNA samples with RIN value ≥ 7
   were used for small RNA enrichment.

Small RNA Library Construction and Sequencing

   Ten small RNA libraries were generated using a TruSeq™ Small RNA Sample
   Prep Kits (Illumina, San Diego, CA, USA) and then sequenced using an
   Illumina HiSeq^TM4000 sequencer (Illumina, San Diego, CA, United
   States) at the Gene Denovo Biotechnology Co., Ltd (Guangzhou, China).
   The clean reads were obtained by removing the reads containing
   adapters, low quality reads with quality scores < Q20 (the proportion
   of read bases whose error rate is <1%) or with unknown nucleotides, and
   the reads shorter than 18nt in length in the raw reads, using fastp
   v0.18.0. First, the clean reads were mapped to GenBank database v209.0
   and Rfam database v11.0 to annotate and remove other non-coding RNAs,
   including ribosomal RNA (rRNA), transfer RNA (tRNA), small nucleolar
   RNA (snoRNA), small nuclear RNA (snRNA), and small cytoplasmic RNA
   (scRNA). Secondly, the clean reads were mapped to the Caprine Genome
   Assembly ARS1
   ([65]ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/704/415/GCF_0017044
   15.1_ARS1) to remove their exons, introns and repeated sequences.
   Subsequently, the remaining clean reads were searched against miRbase
   v22.0 to annotate known caprine miRNAs and known miRNAs from other
   species (named known species-conserved miRNA). Finally, for the reads
   that were not annotated to miRBase V22.0, but matched the Caprine
   Genome Assembly ARS1, they were used to predict novel miRNAs using the
   miReap v.0.2. To ensure the uniquely annotated results for the reads,
   the following annotation ranking was used: rRNA > caprine miRNA >
   caprine miRNA edit > species-conserved miRNA > repeat sequence > exon
   sequence > novel miRNA > intron sequence.

Differentially Expressed miRNAs Analysis and Small RNA Sequencing Results
Validation

   The expression level of miRNAs was first normalized using transcripts
   per million (TPM). The TPM value is calculated by actual reads of each
   miRNA^*10^6 by total reads of all miRNAs. The DESeq v2.0 ([66]22) was
   used to screen differentially expressed miRNAs in Longissimus dorsi
   muscle tissues between LC and ZB goats, using the thresholds of a |fold
   change| > 2.0 and P-value < 0.05. To validate the accuracy of small RNA
   sequencing results, 23 differentially expressed miRNAs were selected
   for reverse transcription-quantitative PCR (RT-qPCR) analysis,
   including eight up-regulated miRNAs (miR-628-5p, miR-885-3p,
   novel-m0312-3p, miR-1994-3p, miR-67-3p, miR-278-3p, miR-307-3p, and
   novel-m0298-5p) and 15 down-regulated miRNAs (miR-381, miR-127-3p,
   miR-200c, miR-136-3p, miR-487b-3p, miR-200a, miR-410-3p, miR-136-5p,
   miR-127-5p, miR-141, miR-200b, miR-276-3p, novel-m0213-5p, miR-2796-3p,
   and miR-429) in Longissimus dorsi muscle of LC goats compared to ZB
   goats. The same RNA samples as those used for the small RNA sequencing
   were used to generate cDNA using a miRNA 1st Strand cDNA Synthesis Kit
   (Accurate Biology, Hunan, China). The caprine U6 and 18SrRNA were used
   as internal references to normalize the relative expression level of