The animal studies were done in accordance with the regulations and guidelines set out by the government of Gansu Province, as well as with the approval of the Animal Health and Ethics Committee of Gansu Agricultural University (Animal Experimentation License No. 2012-2-159).

The experimental animals used in this study have been comprehensively detailed, together with their corresponding management regimens, in previous publications (Zhang et al., 2017b; Zhang et al., 2019). To put it simply, a total of 137 male Hu sheep were obtained from Jinchang Zhongtian Sheep Co., Ltd. (Jinchang China) and transported to Minqin Zhongtian Sheep Co., Ltd. (Minqin China) during the same time frame for the purpose of breeding. The process of weaning was established when the lambs reached 56 days of age. Subsequently, during the initial phases of the experimental study, only lambs displaying resilient development and overall outstanding health were selected as candidates. The lambs were supplied with nourishment in a standardized single pen (0.8 × 1 m), whereby they access to fresh water and food every day until the end of the experiment (180 days of age). The lambs attained an ideal age of 80 days for the start of the official experiment, which was recorded as day one. The performance experiment is concluded when the lambs reach 180 days, hence rendering the official duration of the trial as 100 days. Consequently, all lambs participating in this study had a 2-week transitional period prior before a 10 days pre-feeding phase. Throughout the transitional phase, a regular percentage adjustment was made to the form of feed utilized each day. Additionally, the whole pelleted feed was consumed on the initial day of the pre-fed phase. During the first 10-day pre-fed period and ensuing 100-day official trial, the pellet feed was purchased from Gansu Sanyang Jinyuan Animal Husbandry Co., Ltd., (Gansu China).

Lambs were treated to a feeding regimen structured in 20-day cycles until the completion of the feeding experiment (180 days of age), with their initial weight being measured on the first day of the specified period (80 days of age). The lambs underwent daily weighing before to feeding, using a calibrated electronic scale. No modifications were made to the participants or the equipment utilized over the whole period of the experiment. Furthermore, each sheep was weighed for remaining feed before each weighing period, which was used to calculate feed intake and RFI. For the computational model used in this study, the main reference was the formula of Zhang et al. (2017b). The specific formula used in this particular instance was as follows: Y k = β 0 + β 1 MBW k + β 2 ADG k + e k , MBW = BW i + BW f / 2 0.75 , ADG = BW f ‐ BW i / N , where Y_k represents the average daily feed intake of the ith individual; β₀ regres-sion intercept; β₁ regression coefficient for mid-test metabolic body weight (MBW); β₂ regression coefficient for average daily gain (ADG); e_k represents uncontrolled error of the ith individual; BW_f represents final body weight; BW_i represents initial body weight; and N, experimental period (days).

The methodologies used for the collection and processing methods of the tissues were cited from previous scholarly investigations (Zhang et al., 2017b). The methodology may be concisely summarized in the following manner: all lambs are uniformly transferred to a professional slaughterhouse after the end of the measurement. The RFI values of the six sheep used in this study are provided fully, concerning prior research conducted by our research group. The lambs were subjected to a 24 h of fasting period before to being weighed and then executed in standard procedures. Each liver sample was immediately collected after slaughter process and then preserved in liquid nitrogen for temporary storage. Following the process of slaughter, it was transferred to −80°C for long-term storage until RNA was extracted. As described in previous study, we selected 3 High_RFI and 3 Low_RFI sheep from 137 male Hu lambs for total RNA extraction (Zhang et al., 2017b; Zhang et al., 2019). The total RNA extraction was performed using the TRIzol Reagent (Invitrogen, Waltham, MA, United States) method as per the provided instructions. The NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States) was used to quantify the purity of RNA, while the integrity of RNA was performed using the Agilent Bioanalyzer 2,100 system’s RNA Nano 6000 Assay Kit (Agilent Technologies, CA, United States).

The small RNA library preparation kits were used to produce small RNA sequencing libraries (Galina-Pantoja et al., 2006). There are many steps that involved, as seen below: firstly, whole RNA molecule was used as a template to directly connect the 3′ and 5′ ends of the small RNA with the synthetic adaptors. The synthesis cDNA from total RNA was conducted with M-MuLV reverse transcriptase (NEB, United States), followed by amplification of the resulting cDNA in accordance with the recommended protocols provided by Illumina. The PCR products underwent purification and recovery processes using 8% polyacrylamide gels. These gels exhibited the capability to effectively separate DNA fragments with sizes 140 to 160 base pairs. Following that, the DNA fragments that had undergone purification were dissolved in 8 μL of elution solution. In the end, the evaluation of the library’s integrity on the Agilent Bioanalyzer 2,100 system may be conducted by using DNA high-sensitivity chip. Once the library has undergone qualification, the product was subjected to sequencing on the Illumina HiSeq 2,500 platform, resulting in the generation of a 50 bp single-ended read.

The raw data collected by the sequencing equipment was then given a Base Calling analysis with the intention of producing FASTQ files. The acquisition of clean data included the elimination of extraneous information from the raw data via the use of a customized Perl script. Concurrently, Q20, Q30, and GC-contents of the raw data were obtained. Then, for all subsequent analysis, choose certain range of length from clean reads. The Bowtie (Langmead et al., 2009) method was used for the purpose of aligning small RNA tags with the reference sequence, facilitating the assessment of their expression levels on the reference. Then used miRbase 20.0 as a reference to mapping known miRNAs. In order to exclude protein-coding genes, repetitive sequences, rRNA, tRNA, snRNA, and snoRNA, custom scripts are used to extract miRNAs of predetermined length. Following, the software tools miREvo (Wen et al., 2012) and mirdeep2 (Friedländer et al., 2012) were used to map sequences with the sheep reference genome (Oar_v1.0) to provide predictions about novel miRNAs. The use of the whole rRNA ratio served as an indicator of sample quality. In animal samples, this ratio should ideally be below 40%. Additionally, the cumulative p values for RNA folding were used as a metric for output measure.

The miRNA expression levels were estimated by TPM (transcript per million) (Zhou et al., 2010). The processing procedure differs according on the sample’s quality. To perform differential expression analysis on Low_RFI and High_RFI samples with biological recurrence, use the DESeq R package (version 1.8.3). The filtering conditions for differential expression of miRNAs in this study were p-value <0.05 and |log2 (foldchange)|≤ 0.5 (Tang et al., 2007). The threshold for substantial differential expression is set to the default value. The target gene of miRNA was then predicted for animals using miRanda (Enright et al., 2003). Statistics calculations were carried out to analyze the expression levels of the most prevalent known miRNAs and novel miRNAs in both experimental groups. Additionally, the most common miRNA or DE_miRNA expression was estimated in relation to the anticipated mRNA expression of its target gene.

The target gene candidates of DE_miRNAs were analyzed using Gene Ontology (GO) enrichment analysis (“target gene candidates” in the follows). For GO enrichment analysis, a GOseq-based Wallenius non-central hyper-geometric distribution (Young et al., 2010) was used, which might correct for gene length bias. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) was a database used to understand the advanced functions and benefits of biological systems based on molecular-level information, especially genome sequencing and other highly experimental technologies (http://www.genome.jp/kegg/). To assess the statistical enrichment of target gene candidates in KEGG pathways, we employed the KOBAS (Mao et al., 2005) tool.

Investigate the possible association between DE_miRNA found in this study and Zhang et al. describes DE_mRNA (Cigrovski Berkovic et al., 2020). The construction of a miRNA-mRNA interaction network was facilitated by using Cytoscape software (Shannon et al., 2003). This network was established by including DE_miRNAs and DE_mRNAs based on their specific roles. Specifically, mRNAs exhibiting discernible association with miRNAs were integrated into the miRNA-mRNA interaction network. The same shape with various colors shows types that are up-or downregulated in DE_miRNAs and DE_mRNAs.

3.2 Illumina sequencing generated more than 10 M (million) high quality raw reads for each of the two groups of Hu sheep (except for High_RFI No.1 Hu sheep) (Table 1). The raw data obtained from sequencing is given to a filtering process so as to generate clean reads. Reads with more than 10% N content were removed first. An average of 124 reads was removed in the Low_RFI group, less than 0.01%, and similarly an average of 133 reads was removed in the High_RFI group, less than 0.01% (Supplementary Table S1). Further, after removing low-quality readings, the Low_RFI group averaged 0.32% and the High_RFI group 0.35% being removed (Supplementary Table S1). The number of deletions resulting from the existence of 5 adapter contamination and the presence of ployA/T/G/C was low, with an average 0.00% and 0.03% in the Low_RFI group 0.00% and 0.06% in the High_RFI group, respectively (Supplementary Table S1). More reads were deleted due to the absence of 3 adapter null or insert null, about 0.97% and 0.95%. All samples were filtered to retain clean reads of 98% or more. Finally, the clean reads of each sample were screened for sRNAs within a certain length range (18∼35 nt) for subsequent analysis (Table 1). The length-screened sRNAs were localized to the Ovis aries reference genome to analyze the distribution of small RNAs (Table 1). From the clean data, a total of 10,278,758, 10,014,172, 9,533,396, 8,699,934, 9,391,798, and 13,765,901 mapped reads from the LR1 (Low_RFI No.1 Hu sheep), LR2 (Low_RFI No.2 Hu sheep), LR3 (Low_RFI No.3 Hu sheep), HR1 (High_RFI No.1 Hu sheep), HR2 (High_RFI No.2 Hu sheep), and HR3 (High_RFI No.3 Hu sheep) libraries were retrieved, with over 90% mapping to the Ovis aries reference genome (Table 1). In terms of miRNA expression level, our research results indicate that genes with a TPM <60 retrieved from RNA-seq accounted for about 75% of the total, whereas high-expressed genes, that is, genes with TPM >60, accounted for around 25% (Table 2). However, HR2 accounted for only 7.81%. Screening miRNAs ranged from 18 to 35 nt in length (Figures 1A, B). Furthermore, a variety of non-coding RNAs (ncRNAs) were discovered, including transfer RNAs (tRNAs), snRNAs and miRNAs (Figures 1C, D). Among the identified ncRNAs, a minute fraction constituted recently discovered miRNAs.

We identified 121, 120, and 122 known miRNAs in the High_RFI samples, and 119, 83, and 128 known miRNAs in the Low_RFI samples (Table 3). In all of these reads, approximately 59% of known miRNAs were detected in all samples (Supplementary Table S2). The miRNA the highest abundance across all samples was oar-miR-148a, with mean paired reads from HR1, HR2, HR3, LR1, LR2, and LR3 of 2,885,255, 1,771,013, 2,178,809, 3,049,001, 7,921,051, and 1,323,994, respectively. Among the top 20 expressed miRNAs in each group, the expression levels of 10 miRNAs such as oar-miR-148a, oar-let-7f, oar-miR-143, oar-miR-30a-5p, oar-miR-26a, oar-miR-21, oar-let-7g, oar-let-7i, oar-miR-30d and oar-miR-99a accounted for an average of more than 94% (Supplementary Table S3). The expression analysis shows the top 20 highly expressed miRNAs in the study samples from each respective group (Supplementary Figure S1).

The hairpin structure that is characteristic of miRNA precursors may be used as a means to forecast the existence of novel miRNAs. We identified 73, 90, and 97 novel miRNAs in the High_RFI samples, and 79, 28, and 90 novel miRNAs in the Low_RFI samples (Table 3). Among the novel that were identified, only 23 were identified in all samples (Supplementary Table S3). Of the 117 unique novel miRNAs, novel_31, novel_50 and novel_41were the most expressed in the samples with an average of 1,826, 1,246 and 1,173 reads aligned to these miRNAs, respectively (Supplementary Table S3). All detected new miRNAs were supplied (Supplementary Table S3), and the top 20 expressed novel miRNAs for each group were presented (Supplementary Table S4).

The sequencing data has been submitted to the NCBI Sequence Read Archive (SRA) database under the biological project PRJNA813431. Differential miRNA expression analysis between Low_RFI and High_RFI groups from the same population of sheep with different phenotypic values. The relevant phenotypic data for the selected sheep were detailed in previous studies (Zhang et al., 2022). Phenotypic differences were not significant except for FCR, RFI and feed intake (FI). In total, an average of 11.62 million raw read were obtained from each sample. A total of 247 miRNAs were detected in 6 liver samples, of which four miRNAs (one known miRNA and three novel miRNAs) were identified as differentially expressed (p < 0.05) (Table 4). The Venn diagrams show miRNAs that are uniquely expressed or co-expressed in different groups, with 205 miRNAs co-expressed in both groups (Figure 2A). Four miRNAs were differentially expressed in the Low_RFI group compared to the High_RFI group, including three upregulated and a downregulated (p < 0.05) (Table 4) (Figure 2B). Of all the DE_miRNAs identified, novel_41 and novel_115 were expressed in all samples. The novel_41 was upregulated in Low_RFI group, while novel_115 showed downregulated compared to the High_RFI group (Table 4).

The target genes of ten most highly expressed miRNAs (seven known and three novel) were predicted in the two groups for further functional analysis (Supplementary Table S5). The majority of these target genes mainly involved in various biological processes, including glycolysis, protein synthesis and catabolism, cell growth and proliferation, scavenging of free radicals, as well as cell death and survival (Supplementary Figure S2). In terms of glycolytic processes, the target genes were mainly involved in nucleoside diphosphate kinase activity, nucleoside kinase activity, adenylate kinase activity, hexokinase activity, and other glucose binding activities. For protein synthesis and catabolism, target genes were involved in the positive regulation of protein secretion, protein polymerization and protein deubiquitination, among other roles. For cell growth and proliferation, target genes were involved in platelet alpha granulation, cell cortex proliferation, and the proliferation and development of microtubule cell ribosomes (Supplementary Figure S2).

To further understand the biological functions and roles of these four DE_miRNAs, target genes were identified (miRDB) for highly differentiated miRNAs between Low_RFI and High_RFI gruops. Upregulated miRNAs in Low_RFI group were integrated to several target genes: top ranking genes were DZANK1, CYP26B1, IDH3G, TRAC, IL36RN, UBE2Z, ZMYND12, PDGFD, VSIG2, and ELK4, whereas top-ranking target genes with downregulated miRNAs were CFAP221, HLA-DOB, TOP2B, ACSL4, FRYL, RAB3GAP2, TSPAN9, ILKAP, USP13, and PRR36. The provided diagram illustrates the top 20 anticipated target genes for the four DE_miRNAs (Figure 3). To further elucidate the functions of the DE_miRNAs, we performed enrichment analysis of their candidate target genes. GO enrichment results showed that these target genes were mainly related to metabolism and binding: Biological processes: metabolism, organic metabolism, primary metabolism and macromolecular metabolism. Cellular components: intracellular, intracellular fractions, organelles and membrane-bound organelles. Molecular functions: binding, protein binding, catalytic activity and Hydrolytic enzyme activity (Figure 4A). It was shown that the metabolism in the liver plays an important role in the efficiency of animal feed. KEGG pathway showed DE_miRNAs target genes were mainly enriched in transcriptional dysregulation in herpes simplex virus infection, leishmaniasis, and biosynthesis of amino acids (Figure 4B).

To fully understand the potential RFI effects of miRNA, we used DE_miRNA and their targets genes to create an interactive population network. A total of 1423 DE_miRNAs target genes were identified. We previously discovered 101 DE_mRNAs between the low and high RFI sheep groups using the same objective as the present investigation (Zhang et al., 2022). Among these, seven miRNAs were co-expressed (Figure 5A). However, some target genes were predicted to be the target genes of a single miRNA, and it was observed that upregulation of DE_miRNAs did not always result in downregulation of DE_mRNAs in liver tissue obtained from animals with Low_RFI. Most DE_miRNAs were predicted to primarily target a single differential target gene. As shown in figure, the upregulated miRNA novel_171 targets the upregulated target gene RTP4, and the upregulated miRNA novel_41 targets the downregulated target gene CD274 (Figure 5B). However, the upregulated miRNA oar-miR-485-3p targeted three different genes, the downregulated target gene OAS1 and the two upregulated target genes SHISA3 and PLEKHH2.