Twelve healthy Southdown × Hu F1 sheep were castrated in Qinghuan Mutton Sheep Breeding Company in Gansu Province. The age of the sheep was obtained from sheep breeding records. There were 3 sheep at 0 months old (newborn, namely M0–1, M0–2, M0–3), 3 sheep at 3 months old (sexually immature, namely M3–1, M3–2, M3–3), 3 sheep at 6 months old (sexually mature, namely M6–1, M6–2, M6–3), and 3 sheep at 1 year old (adult, namely Y1–1, Y1–2, Y1–3). We removed testes from 12 sheep after anesthesia and then stored them in an RNA/DNA sample protector (Servicebio, Wuhan, China). The testis of each sheep was dissected longitudinally, and the right testicular tissue of each sheep was collected, part of which was immediately frozen in liquid nitrogen and stored at −80°C for total RNA and protein extraction. The remaining portion will be fixed in a 2.5% glutaraldehyde solution and subsequently processed for paraffin embedding and sectioning.

After total RNA was extracted by Trizol reagent kit (Invitrogen, Carlsbad, CA, United States), the RNA molecules in a size range of 18–30 nt were enriched by polyacrylamide gelelectrophoresis (PAGE). Then the 3′ adapters were added and the 36–48 nt RNAs were enriched. The 5′ adapters were then ligated to the RNAs as well. The ligation products were reverse transcribed by PCR amplification and the 140–160 bp size PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeq Xten by Gene DenovoBiotechnology Co. (Guangzhou, China)

Basic reads are converted to sequence raw data by base calling. Low-quality reads were filtered to remove reads containing 5′ primer contaminants and poly (A). Reads without 3′ adapters and insert tags, as well as Reads shorter than 15 nt or longer than 41 nt in the raw data, were filtered out to obtain clean Reads. For preliminary analysis, the length distribution of clean sequences in the reference genome was determined. Non-coding RNAs are labeled as RNAs, tRNAs, small nuclear RNAs (snRNAs), miRNA, gene, rep, miRNA and unannotation. These RNAs were aligned and subsequently Bowtie (16) searches were performed against Rfam v.10.1¹ (17). Known miRNAs were identified by aligning sequencing reads to the miRBase v22 database² and the sheep reference genome (Qar Rambouillet v1.0). The known miRNA expression patterns in different samples were analyzed (18). Subsequently, the unannotated reads were analyzed by mirdeep2 to predict novel miRNAs (19). The hairpin structure of pre-miRNA and the miRBase database were used to identify the corresponding miRNA star sequence and miRNA mature sequence. The expression of known and novel miRNAs was analyzed using TPM (transcripts per million; miRNAs normalized by TPM) (20). Differential expression analysis of miRNAs between the two groups was performed using the DEG algorithm in the R package (21). The P- values were adjusted by the method of Benjamini and Hoch-berg to control the false discovery rate (FDR). Differentially expressed miRNAs were defined when the adjusted FDR was <0.05. In addition, we used the mean TPM value to calculate fold change (FC) between groups, defining up-regulated (log2FC ≥ 1) and down-regulated (log2FC ≤ 1) miRNAs, respectively. The target of DE miRNA was predicted using Miranda software with the following parameters: single residue pair matching score ≥ 150, ΔG ≤ −30 kcal/mol, and strict 5′ seed pairing required (22). Considering the hypergeometric distribution, R was used for Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DE miRNA-target-get genes (23–25).

Twelve known sheep miRNAs were randomly selected for analysis, and qRT-PCR was used to validate the RNA-seq data. Real-time PCR was performed using a LightCycler 480 II real-time PCR apparatus (Roche, Swiss). Each sample was analyzed in triplicate. microRNA specific primer sequences were designed and synthesized in the laboratory by Qingdao Biotechnology Company based on microRNA sequences obtained from the miRBase database (Release 20.0). Primer information is given in Supplementary Table S1.

We used Illumina Hiseq 2500 sequencing to analyze small RNA populations in 12 libraries obtained from testis of 0-month-old (M0), 3-month-old (M3), 6-month-old (M6), and 1-year-old (Y1) sheep. The results of miRNA sequence quality control are presented in Supplementary Table S2. A total of 157.27 Mb raw reads were obtained. After removing low-quality sequences, aptamers and discarding sequences shorter than 18 nt, 35.58 Mb, 32.67 Mb, 34.01 Mb and 41.01 Mb clean reads were obtained from M0, M3, M6 and Y1 libraries, respectively, for further analysis. 922,960, 1,780,035, 3,372,847, and 4,070,428 unique srRNAs were extracted from M0, M3, M6, and Y1 testes and mapped to the sheep reference genome (Supplementary Table S3). All clean reads were aligned to the miRBase database and recorded as one of the known RNA categories based on their biogenesis and annotation (Figure 1A). As shown in Figure 1A, Y1(7.21%) accounted for the largest proportion of known miRNAs, while M0(6.56%) and M3(5.85%) accounted for the second, and M6(3.3%) accounted for the smallest. However, the highest proportion of unannotated small RNAs was found in the M6(96.13%) and Y1(91.89%) libraries, which represent other types of small RNAs such as piRNAs. A principal component analysis (PCA) of all mapped genes showed that M0, M3, M6, and Y1 group could be distinguished by part along the axis of the first principal component (Figure 1B). Analysis of the read length distribution of all small RNA libraries showed that the dominant length of small RNAs was 22 nt, accounting for at least 37.72% (Figure 1C).

To identify miRNAs in sheep testes, the sequences obtained after removal of small RNAs, including tRNA, snRNA, and rRNA, were compared with the miRBase database to obtain known miRNAs. Sequences that were not annotated into the miRBase database were predicted as novel miRNAs. A total of 787 known miRNAs and 415 novel miRNAs were identified across the 12 libraries (Supplementary Table S4). The frequency of expression of each miRNA in the 12 libraries varied widely, ranging from a few to hundreds of thousands of sequence reads. The expression of these novel miRNAs was very low, ranging from 18 to 25 nt in length with a distribution peak of 22 nt. Overall, the expression of most known miRNAs was highest in the testis of 3-month-old sheep, followed by 6-month-old and finally 0-month-old sheep (Figure 2A). The most highly expressed members in all libraries were members of the mir-199-x, mir-125-x, miR-99-z, miR-151-x, and miR-202-x families, each with more than 100,000 reads.

Using the differential expression genes (DEGs) algorithm from an R package, the transcriptional changes of microRNAs (miRNAs) during sheep testis development were analyzed. Comparisons of miRNA expression profiles between different libraries were presented in Supplementary Tables S5–S10, visualized through scatter plots. Additionally, a scatter plot illustrated the pattern of DE miRNAs among all control groups, revealing that only 138 DE miRNAs were identified between M6 and Y1 testes, albeit with low expression levels (Figure 2B). Between M0 and M3, a total of 217 DE miRNAs were identified, with 139 upregulated and 78 downregulated. Between M0 and M6, 254 DE miRNAs were found, comprising 205 upregulated and 49 downregulated. From M0 to Y1, a remarkable change of 405 DE miRNAs was observed, with 319 upregulated and 86 downregulated. During the transition from M3 to M6, 130 DE miRNAs were identified, of which 97 were upregulated and 33 were downregulated. Lastly, between M3 and Y1, 305 DE miRNAs were detected, with 236 upregulated and 69 downregulated (Figure 2C).

Using the Miranda software, the target genes of miRNAs were predicted. Subsequently, the potential functional roles of these differentially expressed miRNAs were investigated. Potential targets for the differentially expressed miRNAs were screened based on criteria of a minimum single-residue pair total match score of ≥150 and a total energy of ≤ − 30 kcal/mol. Among the 217, 254, 405, 130, 305, and 138 miRNAs identified between M0 vs. M3, M0 vs. M6, M0 vs. Y1, M3 vs. M6, M3 vs. Y1, and M6 vs. Y1, respectively, 5,233, 8,769, 23,588, 696, 8,021, and 180 targets were predicted in 885, 1,679, 2,679, 388, 2006, and 171 target genes, respectively. For the majority of differentially expressed miRNAs, multiple distinct target genes were present; however, for some differentially expressed miRNAs, only a single target gene was identified. Additionally, certain target genes were targeted by multiple differentially expressed miRNAs. For instance, oar-miR-431 was predicted to target 914 genes, miR-199-y was predicted to target 1,510 genes, while BMP4 was uniquely targeted by miR-142-y.

To gain a deeper understanding of the functions of DE miRNAs in sheep testis development, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on the candidate target genes of all DE miRNAs. The GO analysis results are presented in Figure 3. Between M0 and M3, target genes were significantly enriched in categories such as animal organ development, cell morphogenesis, cell part morphogenesis, and Rab GTPase binding (Figure 3A). Between M0 and M6, target genes were significantly enriched in regulation of cellular biosynthetic process, dendrite morphogenesis, tube development, organ morphogenesis, and cell–cell adherens junction (Figure 3B). Between M0 and Y1, target genes were notably enriched in reproductive process, apical junction complex, apical plasma membrane, and ATP binding (Figure 3D). Between M3 and M6, target genes were significantly enriched in cell differentiation, cell morphogenesis involved in differentiation, Wnt signaling pathway, and apical junction complex (Figure 3C). Between M3 and Y1, target genes were enriched in tube morphogenesis, regulation of organelle organization, plasma membrane, Golgi complex, and microtubule binding (Figure 3F). Lastly, between M6 and Y1, target genes were significantly enriched in cell migration, dendrite morphogenesis, organ morphogenesis, and cell–cell junction organization (Figure 3E).

KEGG Pathway Annotation identified significant enrichment of target genes of DE miRNAs in various comparisons. In M0 vs. M3, 30 pathways were enriched, including the Ras, MAPK, Rap1, cAMP, and Notch signaling pathways (Figure 4A). Four key pathways related to testis development were cAMP, MAPK, PI3K-Akt, and Wnt signaling. In M0 vs. M6, target genes were enriched in 30 pathways, such as adherens junction, tight junction, spinocerebellar ataxia, phospholipase D signaling, and endocytosis (Figure 4B). Four testis development-related pathways were MAPK, cAMP, FOXO, and GnRH signaling pathways. Between M3 and M6, 30 pathways were enriched, including Notch, Ras, thyroid hormone, and ErbB signaling (Figure 4C). Five key pathways were Notch, ErbB, endocytosis, Ras, and thyroid hormone signaling.

In M0 vs. Y1, 49 pathways were enriched, including Notch, Ras, thyroid hormone, and ErbB signaling (Figure 4D). Five testis development-related pathways were FOXO, ErbB, endocytosis, Ras, and thyroid hormone signaling. Between M6 and Y1, 30 pathways were enriched, such as Hippo, adherens junction, cAMP, Notch, and Ras signaling (Figure 4E). Five key pathways were Notch, Hippo, cAMP, endocytosis, and FOXO signaling. Similarly, in M3 vs. Y1, 30 pathways were enriched, including tight junction, endocytosis, Ras, adherens junction, and spinocerebellar ataxia (Figure 4F). Four testis development-related pathways were Ras, MAPK, and cAMP signaling (Figure 5).

We performed co-expression cluster analysis of the differentially expressed miRNAs screened from each comparison group and screened the top 5 miRNAs with degree values. After comparison, the number of miRNA target genes in M0 vs. M3 group, M0 vs. Y1 group and M0 vs. M6 group was higher than that in other groups. Notably, oar-miR-29b, oar-miR-143, oar-miR-370-3p, oar-miR-19b, oar-miR-432, and oar-miR-29b were repeated for different times in the 6 comparison groups. TGF family genes related to testicular development were all regulated by differential miRNAs in the five comparison groups (Figure 5).

Randomly selected 12 known sheep miRNAs, including oar-miR-152, oar-miR-10b, oar-miR-21, oar-miR-22-3p, oar-miR-25, oar-miR-99a, oar-miR-379-3p, oar-miR-409-5p, oar-miR-1179, oar-miR-132, oar-miR-374a, and oar-miR-28-y, were used to validate the RNA-seq data. As shown in Figure 6, the qRT-PCR results exhibited similar expression trends compared to the small RNA-seq results.