The Alpine Merino sheep used in this experiment originated from the Gansu Sheep Breeding and Promotion Station. All experimental animals were adult ewes from the same population. According to the difference of wool fib re diameter, coefficient of variation (CV), and mean fibre diameter (MFD): SF group (MFD = 17.68±0.26 μm, CV<20%), EF group (MFD = 19.15±0.37 μm, CV<22%), and M group (MFD = 22.51±0.43 μm, CV<23.6%). Comparisons between the EF and SF (EF–SF), M and EF (M–EF), and M and SF (M–SF) groups were set up by collecting skin tissue from the left shoulders of three sheep from each group. The skin tissue samples (collection area: ~ 2 cm×3 cm) were placed in a freezing tube, immediately frozen in liquid nitrogen, and stored at −80°C.

Total RNA was extracted using the TRIzol reagent (Invitrogen) in accordance with the manufacturer’s protocol. RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent). The libraries were sequenced on an illumina Novaseq 6000 platform.

The fastp tool was used to process raw reads and obtain clean reads for each sample [14]. Which were retained for subsequent analyses. The clean reads were mapped to the reference genome (Oar_v4.0) using HISAT2 [15]. The fragments per kilobase of transcript per million mapped fragments (FPKM) [16] of each gene were calculated, and the read counts for each gene were obtained using HTSeq-count [17]. Differential expression analysis was performed using the DESeq2 tool [18]. A p value<0.05 and a fold change (|FC|)>1.5 were set as the thresholds for identifying significant DEGs. Using the hypergeometric distribution to analyze GO and KEGG significantly enriched items, and screen for a threshold of p<0.05.

Small RNA libraries were constructed, purified, and sequenced using the NEBNext Small RNA Library Prep Set for Illumina kit (NEB#E7330S; NEB). Filtering of low quality reads, The length distribution of the clean sequences in the reference genome was determined, then the sequences were aligned and subjected to the Bowtie [19] search against Rfam v.10.1, ribosomal RNA, small cytoplasmic RNA, cis-regulatory RNA, small nuclear RNA, transfer RNA, and other RNAs were annotated and filtered. Next, cDNA sequences and species-specific repeat sequences from the Repbase database [20] were identified with Bowtie software. Mature miRNAs were identified by aligning against the miRBase v22 database, and the expression patterns in different samples were analyzed. Unannotated reads were analyzed by miRDeep2 [21] to predict novel miRNAs. Based on the hairpin structure of each pre-miRNA and the miRBase database, the corresponding miRNA star and miRNA mature sequence were also identified. Differentially expressed miRNAs (DEmiRNAs) were calculated and filtered with the threshold of p<0.05 and |FC|>1.5. The targets of DEmiRNAs in animals were predicted using miRanda software [22], with the following parameter score (S)≥150, a change in free energy (ΔG)≤−30 kcal/mol, and strict 5′ seed pairing requirements. GO enrichment and KEGG pathway enrichment analysis of DEmiRNAs target Genes were performed using R, based on the hypergeometric distribution.

StringTie was used to count reads for determining the original expression levels of lncRNA, which were normalized using FPKM. We screened for lncRNAs long than 200 bp and coverage exceeding 3. the selected lncRNAs were using PLEK, CNCI, and Pfamscan software to identify those with high credibility. To identify circRNA, we utilized find_circ. The 20 bp sequences at both ends of the unaligned reads from the HISAT2 alignment results were extracted and employed as anchor sequences. These sequences were realigned to the genome using Bowtie2 for circRNA detection. The read count values obtained from the circRNA identification results of find_circ were used as the original expression level of circRNAs, which were normalized using transcripts per million.

DESeq was used to normalize lncRNAs and circRNAs and to calculate FCs, distinguishing differentially expressed lncRNAs (DELs) from differentially expressed circRNAs (DECs) based on the criteria of |FC|>1.5 and p-value<0.05. The k-means clustering algorithm was employed to analyze the DELs and DECs, which were functionally analyzed using the GO and KEGG databases.

Skin samples from different groups were ground to a powder in lysis solution (Biyuntian) containing phosphatase inhibitor (Roche) and protease inhibitor PMSF (Biyuntian) The samples were then centrifuged to determine the protein concentration and molecular mass in the supernatant, and adjusted to the same concentration and volume. The samples were incubated in DTT, to which an appropriate volume of iodoacetamide was added and mixed well. Acetone was added to precipitate the protein, and the mixture was centrifuged to collect the precipitate, which was redissolved in NH₄HCO₃ (Source Leaf Bio). Trypsin Trypsin-TPCK (Wallis) was added at a concentration of 1/50th of the sample mass, incubated overnight at 37°C for digestion, and terminated by adjusting the pH to ~3 with phosphoric acid. The digested peptides were desalted using SOLA SPE 96-well plates, separated on an 1100 HPLC System (Agilent) using a Nano Chrom-C18 column (5 μm, 150 mm×2.1 mm) and lyophilized for MS. Peptides were separated in 90 min at a flow rate of 300 nl/min on a 25 cm×75 m column (1.6 μm C18, ionopticks). Mobile phases A and B were 0.1 vol% formic acid solution and 80:20:0.1 vol% ACN: water: formic acid, respectively. The total run was 60 min (45 min, 5%–27% B; 45–50 min, 27%–46% B; 50–55 min, 46%–100% B; 55–60 min, 100% B). Liquid chromatography was coupled online to a hybrid TIMS quadrupole TOF mass spectrometer (Bruker timsTOF Pro) via a Captive Spray nano-electrospray ion source. The conditions of DDA mass spectrometry are shown in the attached Materials and Methods Supplement 1. To perform data independent acquisition (DIA), we used the instrument control software (Bruker otofControl v6) to define quadrupole isolation windows as a function of the TIMS scan time (diaPASEF). The conditions of DIA mass spectrometry are shown in the attached Materials and Methods Supplement 2. The mass and concentration of proteins are shown in Supplement 3.

Spectronaut (ver. 15.3.210906.50606) was used to search all raw data against the sample protein database, performing the database search with trypsin digestion specificity and considering alkylation on cysteine as a fixed modification. The false discovery rates (FDRs) for proteins, peptides, and peptide-spectrum matches were all set to 0.01. For DIA data, the quantification FDR was also set to 0.05, and the quantification was performed at the MS2 level. Proteins with a p-value<0.05 in the three group comparisons (EF–SF, M–EF, M–SF) were considered differentially expressed proteins (DEPs), and were subjected to GO and KEGG analyses.

The cor function in R was used to analyze correlation coefficients (CCs) and p-values between the target gene, DELs, DECs, and DEPs, employing screening conditions of |CC|≥0.9 and p<0.01. The igraph package was employed to construct the ceRNA competitive regulatory network.

Seven RNAs were randomly selected for real time polymerase chain reaction (RT-PCR) analysis, with the ARPC5 gene as an internal reference. The reactions were performed using the PerfectStart Green qPCR SuperMix kit on the LightCycler 480 II Real-Time PCR System (Roche). The expression levels of the genes were analyzed using the 2^−ΔΔCt method. The relationship between the expression levels of these seven genes in different groups and the sequencing data was analyzed using the cor function in R language.

A total of 18,931 mRNAs were identified in this study (Supplement 4), Counting of the miRNAs in the group samples showed that most were between 20–23 bp, with similar expression levels in each sample (Figure 1A). Differential analysis identified 16 DEmiRNAs in the EF–SF comparison, of which eight were significantly upregulated and eight were significantly downregulated, and 16 DEmiRNAs in the M–EF comparison, of which 11 were upregulated. The highest number of DEmiRNAs (24) was identified in the M–SF comparison, of which 15 were upregulated. Only a few of these DEmiRNAs have been previously assigned official names. The prediction of target genes and binding free energy of these DEmiRNAs revealed that target genes could be paired with the seed region sequence of the miRNAs (Figures 1B–1D). Because most of the DEmiRNAs did not have official names, we analyzed their functions on the basis of genomic location. In the EF–SF comparison, the novel76mature miRNA and its target genes were mainly enriched in the PI3K Akt, p53, and Wnt signaling pathways (Figure 1E). In the M–EF comparison, several DEmiRNAs and their target genes were enriched in valuable pathways. Among them, novel47mature and its target genes were found to be involved in the nuclear factor (NF)-kappa B, transforming growth factor (TGF)-beta, and Wnt signaling pathways. Novel76mature was also identified in the M–EF comparison, where it participated in different regulatory pathways (Figure 1F). In the M–SF comparison, six DEmiRNAs were enriched in pathways related to hair follicle development. Multiple target genes of novel333mature showed involvement in multiple pathways, namely negative regulation of the Wnt signaling pathway, the Wnt signaling pathway, and skin development. Other miRNAs, such as novel112mature, novel150star, and novel181_mature and their target genes were only enriched in one hair follicle-related pathway (Figure 1G).

Sequencing data revealed 373,461 and 377 DELs in the EF–SF, M–EF, and M–SF comparisons, respectively. There were large numbers of DECs in all three pairings, with 1,233 DECs in the EF–SF comparison and 1,111 DECs in the M–EF comparison. Additionally, the total number of DECs in the three comparisons was higher than that in the upregulated group (Supplement 5). Next, we applied clustering algorithms to partition DECs and DELs within each group comparison into two distinct clusters, enabling the investigation of expression patterns across experimental groups. In the EF–SF comparison, Cluster 1 DECs were predominantly upregulated in the EF group, while Cluster 2 DECs were predominantly upregulated in the SF group. Functional enrichment analysis revealed that Cluster 1 DECs were significantly enriched in pathways associated with skin development (e.g., epidermal differentiation and hair follicle morphogenesis), with pathway-related DECs demonstrating consistently higher expression levels in the EF group than in the SF group (Figure 2). In contrast, DECs in Cluster 2 were enriched in a broader spectrum of biological processes, with key DECs involved in representative pathways (e.g., extracellular matrix organization and immune response) showing elevated expression in the SF group. For DELs in the EF–SF comparison, Clusters 1 and 2 displayed divergent expression profiles and functional annotations. DELs in Cluster 1 were predominantly associated with the regulation of genetic material synthesis (e.g., DNA replication and RNA processing), while DELs in Cluster 2 were enriched in metabolic processes and cellular communication pathways, including signal transduction and intercellular transport (Supplement 6). In the M–EF comparison, the numbers of DECs and DELs in Clusters 1 and 2 were comparable. DECs in Cluster 1 showed significant enrichment in pathways associated with skin and hair follicle development (e.g., epidermal morphogenesis and Wnt signaling). Notably, DECs within these pathways exhibited consistently higher expression levels in the EF group than in the M group (Figure 3). In contrast, DECs in Cluster 2 were enriched in biologically relevant pathways; however, pathways enriched by DECs in the M-EF comparison group lacked directional specificity and exhibited a broad distribution. For DELs in the M–EF comparison, Clusters 1 and 2 displayed divergent expression patterns between the two groups. Functionally, DELs in Cluster 1 were predominantly involved in biological binding processes (e.g., protein–DNA interaction), while DELs in Cluster 2 participated in both binding processes and metabolic regulation (e.g., lipid catabolism) (Supplement 7). In the M–SF comparison, DECs in Cluster 1 were enriched in pathways involved in TGF-beta signaling and skin development, with pathway-related DECs demonstrating overall higher expression in the M group than in the SF group. DECs in Cluster 2 were enriched in more functionally significant pathways, including hair cycle and hair follicle development, with pathway-associated DECs exhibiting elevated expression in the SF group (Figure 4). DELs in the M–SF comparison retained consistent expression patterns with prior groups, though neither cluster showed a distinct functional orientation (Supplement 8).

Next, we compiled the target genes enriched in the DECs, DELs, and DEmiRNAs associated with the skin follicle development-related pathways in each group to construct a ceRNA regulatory network based on the expression relationships among these RNA types. The EF–SF comparison was mainly dominated by the target genes TP53I3 and PPP2R3A. Furthermore, TP53I3 was negatively correlated with circ_4794, circ_ 0690, TCONS_00009822, and TCONS_00056385, suggestive of competitive relationships (Figure 5A). The target genes in the M–SF comparison were negatively correlated with most of the DELs and DECs, suggesting that multiple circRNAs or lncRNAs may compete for binding to the same miRNA (Figure 5B). DELs were positively correlated with most of the mRNAs in the M–EF comparison; however, PRKCA was negatively correlated with both DECs and DELs, suggesting that these DECs and DELs may share an miRNA with PRKCA (Figure 5C).

Correlation analysis of the transcriptomes and proteomes in the EF–SF, M–EF, and M–SF group comparisons revealed that the CC between mRNA and protein expression in M-SF pairing was the highest (0.41) (Figure 6A), We selected genes and proteins enriched in wool fiber diameter-related regulatory pathways and calculated their CCs. The results showed that the differences in specific DEPs were greatest between the M and SF groups, with very obvious differences in DEGs observed in all three group comparisons (Figure 6B). Analysis of the target genes of skin development-related miRNAs in association with specific DEGs showed that LOC10112037 and LOC105614079 had positive correlations with most of the target genes, suggesting that they may act as key genes to regulate skin development (Figure 6C). Therefore, we constructed the main regulatory network comprising group-specific DEGs, miRNAs, circRNAs, lncRNAs, proteins, and target genes. Key nodes in this network were LOC10112037, LOC105614079, IGFBP1, and IGFBP4. The majority of RNAs exhibited positive correlations with one other, while some target genes, circRNAs, and lncRNAs showed negative correlations, suggesting a complex regulatory relationship among these RNAs (Figure 6D).

We selected seven RNAs, namely oar-miR-23a, oar-miR-150, HOXC13, COL1A1, MAPK9, CHD8 and IGFBP4, and analyzed their expression levels across different groups with varying fiber diameters. The expression levels measured by qRT-PCR showed an extremely strong correlation with the RNA-Seq data (Figure 6), thereby confirming the authenticity and reliability of the RNA-Seq data in this study (Figure 7).