Integrated Analysis of miRNA and mRNA Expression Profiles Reveals Functional miRNA-Targets in Development Testes of Small Tail Han Sheep

Small Tail Han Sheep is a highly valued local breed in China because of their precocity, perennial estrus, and high fecundity. The average annual lambing rate of ewes is as high as 180–270%, the semen of ram has characteristics of high yield, high density, and good motility. To reveal the key miRNAs and miRNA-targets underlying testis development and spermatogenesis in male Small Tail Han Sheep, integrated analysis of miRNA and mRNA expression profiles in 2-, 6-, and 12-month-old testes was performed by RNA-seq technology and bioinformatics methods. The results showed that total of 153 known sheep miRNAs and 2712 novel miRNAs were obtained in 2-,6 - and 12-month-old Small Tail Han Sheep testes; 5, 1, and 4 differentially expressed (DE) known sheep miRNAs, and 132, 105, and 24 DE novel miRNAs were identified in 2- vs. 6-, 6- vs. 12-, and 2- vs. 12-month-old testes, respectively. We combined miRNA results of this study and the mRNA results obtained in our previous study to predict the target mRNAs of DE known sheep miRNAs; 131, 10, and 15 target mRNAs of DE known sheep miRNAs and 76, 1, and 11 DE miRNA–targets were identified in the three groups, respectively. GO and KEGG analyses showed that: in 2- vs. 6-month-olds, the target genes of DE known sheep miRNAs were involved in 100 biological processes and 11 signaling pathways; in 6- vs. 12-month-olds, the target genes of DE known sheep miRNAs were involved in 4 biological processes; and in 2- vs. 12-month-olds, the target genes of DE known sheep miRNAs were involved in 17 biological processes and 4 signaling pathways. Three miR–target regulatory networks were constructed based on these DE miRNA–targets. The key miRNA-Targets involved in testis development and spermatogenesis were screened. 6 known sheep miRNAs and 6 novel miRNAs were selected to validate the accuracy of miRNA sequencing data by qRT-PCR. The binding sites of oar-miR-379-5p with WNT8A was validated by a dual luciferase reporter gene detection system.

sperm density, and good sperm motility. Male Small Tail Han Sheep can be used as a good research subject. Sperm quality is closely related to testis tissue structure, testis development and spermatogenesis are obvious factors that affect male reproductive efficiency (Gong et al. 2013;Zhu et al. 2016). Therefore, understanding the molecular regulatory mechanisms and finding out the molecular marker underlying testis development and spermatogenesis in Small Tail Han Sheep will provide theoretical basis for improving male sheep reproductive efficiency.
With the wide application of molecular biology technology in animal genetics and breeding, it is possible to systematically clarify the molecular regulatory mechanisms underlying gene functions (Miao and Qin 2015). Gene regulation is at the core of all biological behaviors and phenotypes, which is a multi-mode spatial and temporal dynamic process and mainly regulated at transcriptional and post-transcriptional levels . Generally, transcriptome refers to the sum of all RNAs transcribed by an organism in a particular environment or physiological state, including mRNA and noncoding RNA (such as rRNA, tRNA, microRNA, siRNA, piRNA, circle RNA, and lncRNA). miRNAs are important noncoding RNAs that regulate gene expression and participate in many biological and metabolic processes, such as growth, development, reproduction, and disease, by specifically binding to 39untranslated regions (UTRs) of mRNA to degrade a target gene or inhibit target gene translation (Ambros 2004;Ran et al. 2014). Increasing evidence has shown that miRNAs play key roles in testis development and spermatogenesis (Kotaja 2014;Pratt and Calcatera 2016). Testis development and spermatogenesis have obvious stage specificity, miRNA expression changes with testis development and spermatogenic cell development (Zhang et al. 2010;Lian et al. 2012;Wu et al. 2014;Kasimanickam and Kasimanickam 2015;Kasimanickam and Kasimanickam 2015;Tariq et al. 2016). To date, many miRNAs have been identified and studied in different species, such as humans, mice, pig, cattle and chicken (Yan et al. 2007;Zhang et al. 2010;Huang et al. 2011;Luo et al. 2015;Wang and Xu 2015;Wu et al. 2017b). However, miRNAs identification and molecular mechanism exploration has not been investigated in sheep testis as in other livestock species. Therefore, it is necessary to investigate miRNA and mRNA expression profiles to understand the molecular regulatory mechanism and identify the key miRNA-targets underlying testes development and spermatogenesis in Small Tail Han Sheep.
In our previous publication, histological structure and mRNA expression profiles of 2-, 6-, and 12-month-old Small Tail Han Sheep testes were investigated using hematoxylin-eosin staining and RNA-seq technology, respectively (Bai et al. 2017). In this study, we analyzed miRNA expression profiles in 2-, 6-, and 12-month-old Small Tail Han Sheep testes, integrated the miRNA expression profile in this study and the mRNA expression profile obtained in our previous study to construct miRNA-target regulatory networks, selected the key miRNAs and their target mRNAs involved in testis development and spermatogenesis in Small Tail Han Sheep. These results of our study can help to identify the molecular markers affecting reproductive efficiency of male Small Tail Han Sheep and provide theory evidence for further finding out the molecular markers affecting male sheep reproductive ability.

Sample preparation
Two-, 6-, and 12-month-old male Small Tail Han Sheep from the nucleus herd bred by Songyuan Sheng Hua Animal Husbandry Co. Ltd (Jilin, China) were used for this experiment. Three different rams of each age were selected for sample collection. Left testis of each ram was collected by surgical operation and immediately frozen in liquid nitrogen. Intratesticular tissues were used for total RNA extraction. Samples were the same as those used in our previous study (Bai et al. 2017).
miRNA library construction and sequencing Total RNA was extracted from testis tissues using Trizol Reagent (Invitrogen, USA). RNA integrity number (RIN) and RNA quality were determined using Agilent Tapestation 2200 (Agilent, USA) and NanoDrop 2000 (Thermo Scientific, USA) instruments. RNA samples with RIN $6.0 were considered acceptable for miRNA library construction. The miRNA library of each sample was prepared using an Ion Total RNA-seq Kit v2 (Life Technologies, USA) following the manufacturer's instructions. Briefly, 13-34nt RNAs were captured and enriched using magnetic beads with oligo (dt); subsequently, 13-34nt RNAs were fractionated into short fragments using RNase III and an Ion adaptor. The short RNA fragments were ligated with 39 and 59 adaptors. Then, the short mRNA fragments were reverse-transcribed and amplified to synthetize double-stranded cDNA. Emulsion PCR was performed using the miRNA cDNA library as the template. RNA-seq was conducted on an ABI Ion Proton instrument by NovelBio Bio-Pharm Technology Co. Ltd (Shanghai, China). Raw sequencing data were assessed by Fast-QC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc) (Ramayo-Caldas et al. 2012).

DE miRNA analysis
The raw sequencing reads were filtered, and the clean reads were mapped to miRBase 21.0 (http://www.mirbase.org/) to identify known sheep miRNAs (oar-miRNAs). Then, the remaining clean reads that were not mapped to sheep miRBase were mapped to sheep genome Oar3.1 to predict sheep novel miRNAs by miRDeep 2.0.0.8. The remaining clean reads that did not map to sheep genome Oar3.1 were mapped to human, rat, and mouse miRBase using the BWA algorithm to predict novel miRNAs. Only the mapped reads were used for miRNA expression analysis. miRNA expression levels were quantified by the reads per kilobase of transcript per million mapped reads (RPKM) method. The RPKM value was used to calculate expression, and the upper-quartile algorithm was used to correct expression to help produce accurate results for some miRNAs that occurred in low abundance. The DE miRNAs were identified by the EB-Seq algorithm (criteria:|log 2 FoldChange |.0.585; false discovery rate, FDR , 0.05); FDR was calculated to correct the P-value (Wright and Simon 2003).
Venn, Series-cluster and hierarchical clustering analysis of DE miRNAs Venn analysis was performed to visualize the DE mRNA and miRNA changes by stage using Microsoft Excel 2016. Series-cluster analysis was performed to identify the expression profiles of miRNAs using the STEM program (http://www.cs.cmu.edu/$jernst/st/). The DE miRNAs in 2vs. 6-, 6vs. 12-, and 2vs. 12-month-old testes were categorized into 8 expression profiles. Hierarchical cluster analysis was performed with Cluster 3.0 and Tree View 1.6 programs (http://rana.lbl.gov/eisen).

Target prediction of DE miRNAs and integrative analysis
The TargetScan and miRanda algorithms were used to predict putative target mRNAs of DE miRNAs. The putative target mRNAs of DE miRNAs were crossed with DE mRNAs in 2vs. 6-, 6vs. 12-, and 2vs. 12-month-old testes, respectively, which were obtained in our previous study (Bai et al. 2017). Negative correlation analysis was performed to determine the relation of DE miRNAs and mRNAs. Pearson correlation coefficients were used to determine whether the expression levels of each particular miRNA and its target mRNAs were negatively correlated (correlation , 0, P , 0.05, FDR , 0.05). If a miRNA and its target mRNA were negatively correlated, the mRNA was identified as a candidate target of miRNA.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses
The GO database (http://geneontology.org/) was used to analyze the main functions of candidate target mRNAs. Biological process (BP) terms were annotated for candidate target mRNAs. Similarly, the KEGG database (http://www.genome.jp/kegg/) was used to analyze the pathways in which the candidate target mRNAs are involved. Fisher's exact and x 2 tests were used to classify the significant GO terms and pathways of DE mRNAs (P , 0.05), and FDR was calculated to correct the P-value ). All of the candidate target mRNAs were analyzed using the GO and KEGG pathway databases.

miRNA-target regulatory network construction
To identify all possible miRNA-mRNA interactions, miRNA-target regulatory networks were built with Cytoscape 2.8.3. The expression values for DE miRNAs and their candidate target mRNAs were imported into Cytoscape to generate and visualize functional networks.

RNA-seq data validation
Six known sheep miRNAs and 6 novel mRNAs were randomly selected and analyzed to validate RNA-seq data using qRT-PCR. Synthesized cDNA was used as the template for qRT-PCR. qRT-PCR amplification for miRNAs was performed using the miRcute Plus miRNA qPCR Detection Kit (FP411, TIANGEN BIOTECH, China). The 20-ml reaction system contained 10ml 2·miRcute Plus miRNA Premix, 0.4 ml universal reverse primer (200nM) provided by the kit, 0.4 ml miRNA sequence-specific forward primer (200nM), 2ml cDNA, and 7.2 ml ddH 2 O. The reaction conditions were 95°for 15 min, and 40 cycles of 94°for 20s, 60°for 34s; the reaction was completed with a melting curve analysis according to the manufacturer's instructions. qRT-PCR was performed on a Bio-Rad CFX96 system (Bio-Rad, USA). All qRT-PCRs were completed with a melting curve analysis. The 2 2DDCt method and internal control gene U6 was used to calculate miRNA expression levels. Specific primers were designed and synthesized by Sangon Biotech (Shanghai, China); primer information is shown in Table 1.

Plasmid vector construction and transfection
The sequence fragment of WNT8A-39UTR that binds to oar-miR-379-5p (wild-type), WNT8A-39UTR with an oar-miR-379-5p binding site mutation (mutant-type), and positive control (miRNA inhibitor) were designed by Target Scan 7.1 and RNA22 v2.0, and synthesized by total gene synthesis. Sequence information of wild-and mutant-type target gene fragments and the miRNA inhibitor fragment are shown in Table  2. The pmirGLO Dual-Luciferase miRNA Target Expression Vector was used to confirm the putative binding site of WNT8A-39UTR with oar-miR-379-5p. pmirGLO was digested by restriction endonuclease SacI and XhoI. pmirGLO and target gene fragments were ligated by T4 DNA ligase. 293T cells were washed and incubated with 4ml serumfree medium for 24h. The miRNA (or miRNA inhibitor) was mixed with 250ml serum-free medium and 1.6mg recombinant vector, and Lipofectamine 2000 transfection reagent (Life Technologies, USA) was mixed with 250ml serum-free medium. Then, the 2 mixtures were combined and incubated for 25min at room temperature. The Lipofectamine 2000-miRNA mixture was added to the cells and incubated at 37°for 5h. Subsequently, 500mL DMEM containing 10% fetal bovine serum was added to the cells and incubated in 5% CO 2 at 37°. The cells were collected after culturing for 24h and 48h. Then, 7 experimental groups were set up: WNT8A-WT+mimics NC, WNT8A-WT+oar-miR-379-5p, WNT8A-MUT+mimics NC, WNT8A-MUT+oar-miR-379-5p, PC +mimics NC, PC+ oar-miR-379-5p, and a blank group. There were 3 replicates of each group.
Dual luciferase reporter detection Luciferase and Renilla activities were measured using a Dual-Luciferase Reporter Assay Kit (Promega Corporation, USA) after transfection for 24h and 48h, respectively.

Statistical analysis
The data were analyzed using SPSS 18.0, and all results were presented as means 6 SD. The significance of differences between groups were assessed by ANOVA. Differences were considered statistically significant when P , 0.05.
n Table 2 Sequence information of wild-and mutant-type target gene fragments and positive control fragments In Venn analysis, DE mRNA and miRNA changes can be visualized by stage (Figure 1 a b). In series-cluster analysis, we categorized the 955 DE mRNAs and 211 DE miRNAs into 8 possible expression profiles to observe the overall expression tendencies (Figure 1c d). Subsequently, based on these DE mRNAs and DE miRNAs, heat maps were generated by hierarchical clustering analysis (Supplemental figure 1 a b).
Candidate target mRNA GO analysis The GO analysis results are shown in Table S7. In 2vs. 6-month-old testes, the candidate target mRNAs were significantly enriched in 100 BP terms, such as apoptotic process, meiosis I, regulation of acrosome reaction, protein phosphorylation, and epithelial cell maturation; the top 15 BP terms are shown in Figure 2 a: peptidyl-serine phosphorylation, positive regulation of neuron apoptotic process, ionotropic glutamate receptor signaling pathway, epithelial tube morphogenesis, regionalization, pteridine-containing compound metabolic process, regulation of epithelial cell differentiation involved in kidney development, mammary gland development, regulation of transcription involved in anterior/posterior axis specification, polarity specification of Figure 2 Top 15 BP terms of candidate target mRNAs. a. Two-vs. 6-month-old testes; b. 6-vs. 12-month-old testes; c. 2-vs. 12-month-old testes; P , 0.05.

Figure 3
Top 15 pathways of candidate target mRNAs. a. Two-vs. 6-month-old testes; b. 2-vs. 12-month-old testes; red terms represent P , 0.05. proximal/distal axis, canonical WNT signaling pathway involved in neural crest cell differentiation, regulation of renal sodium excretion, positive regulation of apoptotic process by virus, iron ion import, and establishment of organ orientation. In 6vs. 12-month-old testes, the candidate target mRNAs were significantly enriched in 4 BP terms, as shown in Figure 2 b: multicellular organism reproduction, binding of sperm to zona pellucida, single fertilization, and proteolysis. In 2vs. 12-month-old sheep testes, the candidate target mRNAs were significantly enriched in 17 BP terms, such as adenosine metabolic process, nucleotide metabolic process, purine nucleobase metabolic process, cellular response to oxidative stress, and RNA processing; the top 15 BP terms are shown in Figure 2 c: positive regulation of adenosine receptor signaling pathway, negative regulation of caveolin-mediated endocytosis, thiamine metabolic process, negative regulation of pinocytosis, adenosine metabolic process, positive regulation of cell projection organization, regulation of clathrin-mediated endocytosis, regulation of Cdc42 GTPase activity, skeletal muscle fiber development, regulation of sensory perception of pain, nucleotide metabolic process, purine nucleobase metabolic process, regulation of GTPase activity, cellular response to oxidative stress, and regulation of Rho protein signal transduction.
oar-miR-379-5p targets WNT8A According to GO and pathway analysis, 78 target mRNAs of 10 DE known sheep miRNAs in 3 groups were further studied to identify the testis development-and spermatogenesis-related biological processes and signaling pathways in which these target mRNAs were involved. It was found that WNT8A was closely related to tissue development processes, such as multicellular organismal development, cell morphogenesis, and cell fate commitment, and signaling pathways such as WNT, melanogenesis, and Hedgehog signaling pathways; the GO and pathway terms associated with WNT8A are shown in Table S9. Thus, oar-miR-379-5p and WNT8A were selected to verify their targeted regulatory relationship. The binding site of WNT8A-39UTR with oar-miR-379-5p is shown in Figure 6 a. pmirGLO plasmid vector and its working principle are shown in Figure 6 b c. Luciferase assays were performed to confirm the potential binding between oar-miR-379-5p and WNT8A to further validate the interaction of miRNA-mRNA. The luciferase activity of the WNT8A receptor decreased 74.2% and 61.8% upon co-transfection with oar-miR-379-5p mimics (P , 0.05) compared with co-transfection with negative control mimics at 24h and 48h, respectively (Figure 6 d e), which indicated that oar-miR-379-5p can directly target WNT8A-39UTR.

DISCUSSION
Sheep reproductive efficiency is important for livestock husbandry; testis development and spermatogenesis are key factors influencing male sheep reproductive efficiency, many functional genes regulate testis development and spermatogenesis (Gong et al. 2013;Ding et al. 2016). The aim of this study was to identify the key miRNAs and functional miRNA-targets in developing testis of Small Tail Han Sheep using integrated analysis of miRNA and mRNA expression profiles. The analysis will help us understanding the molecular regulating mechanism in testis development and reveal the molecular markers related to reproductive efficiency of male Small Tail Han Sheep. Figure 5 qRT-PCR validation of RNA-seq data for miRNAs. a. RNA-seq data. Values were calculated and normalized by EB-Seq algorithm if the fold changes were .1.5 and FDR , 0.05. b. RT-PCR analysis of 6 novel miRNAs and 6 known miRNAs. The fold change in each time point was relative to 2 months. The expression values were calculated by the 2 2DDCt method. Different capital letters represent extremely significant differences (P , 0.01); different letters indicate significant differences (P , 0.05), whereas the same letters indicate no significant differences (P . 0.05).
In this study, 2-, 6-, and 12-month-old corresponding to infant, adolescent and adult stages in male Small Tail Han Sheep presenting the dynamic process of testis development. Only 5, 1, and 4 DE known sheep miRNAs were identified in 2vs. 6-, 6-vs 12-, and 2-vs 12-month-old testes, respectively, which is because of the limited information in sheep miRBase. Novel miRNAs were predicted from sheep genome, human, rat, and mouse miRBase; therefore, more DE novel miRNAs were identified than DE known sheep miRNAs. DE novel miRNAs can help to enrich our knowledge about the key miRNAs involved in sheep testis development and provide new references for the sheep miRBase. The DE known sheep miRNAs were reported in previous studies, miR -409-3p, miR-665-3p, miR-379-5p, miR-495-3p, miR-23b, and miR-758-3p are all involved in cell proliferation, differentiation, migration, and apoptosis (Weng et al. 2012;Dong et al. 2016;Chen et al. 2017a;Chen et al. 2017b;Jiang et al. 2017;Xie et al. 2017). Based on our results, we can infer that oar-miR-409-3p, oar-miR-665-3p, oar-miR-379-5p, oar-miR-495-3p, oar-miR-23b, and oar-miR-758-3p may play key roles in growth, proliferation, and apoptosis of testis cells. oar-miR-200a, oar-miR-200b, and oar-miR-200c belong to the miR-200 family, which plays a key role in epithelial-mesenchymal transition (Mongroo and Rustgi 2010;Pan et al. 2017). It can be conferred that oar-miR-200a, miR-200b, and miR-200c may be associated with seminiferous epithelium development. Currently, these miRNAs are mainly investigated in the context of invasion and migration of cancer, the specific functions of these miRNAs in sheep testes need to be further explored and validated.
GO and KEGG pathway analyses can help elucidate the functions and pathways that candidate target genes involved in. As our results shown: in 2vs. 6-month-old testes, most of the GO terms were related to spermatogenesis (such as apoptotic process, meiosis I, and acrosome reaction), and most of the pathways were related to tissue development (such as Hedgehog and WNT signaling pathways); in 6vs. 12-monthold testes, GO terms were related to male reproductive processes (multicellular organism reproduction, binding of sperm to zona pellucida, single fertilization, and proteolysis). In 2vs. 12-month-old sheep testes, most GO terms and pathways were related to testis metabolic and biosynthesis processes.
Integrative analysis of miRNA and mRNA expression profiles can help elucidate molecular regulatory mechanisms (Miao and Qin 2015). This study is the first attempt to identify key miRNA-targets by integrated analysis of miRNA and mRNA expression profiles in Small Tail Han Sheep testes. Regulatory networks of DE known sheep miRNAs and their target mRNAs can help us more precisely identify the functions of these miRNAs in testis development and spermatogenesis in sheep. In this study, 3 miRNA-target networks were constructed and several miRNA-targets were identified. According to GO and pathway analysis, WNT8A (the target of oar-miR-379-5p) was determined to be involved in biological processes such as multicellular organismal development, endoderm development, and palate development, and signaling pathways such as the WNT, Hedgehog, and Hippo signaling pathways. It was previously reported that miR-379 is a key miRNA that inhibits proliferation, invasion, and migration of tumor Figure 6 Detection of interactions between oar-miR-379-5p and WNT8A by dual luciferase reporter system. a. Binding site sequence of miRNA and target gene WNT8A. b. Structure of pmirGLO plasmid vector. c. Working principle of pmirGLO plasmid vector. The dual luciferase reporter gene was inserted the 39UTR sequence of the target gene, which can bind with miRNA into the downstream of the luciferase gene in pmirGLO plasmid vector, the translation of luciferase from firefly will be inhibited and the luciferase activity will eventually decrease when endogenous miRNA or introduced exogenous miRNA are combined with the inserted 39UTR sequence. d. Relative luciferase activity after transfecting oar-mir-3799-5p and WNT8A recombinant plasmid for 24h. e. Relative luciferase activity after transfecting oar-mir-3799-5p and WNT8A recombinant plasmid for 48h. cells, and plays an important role in tumor occurrence and development; miR-379 can also inhibit proliferation, invasion, and migration of vascular smooth muscle cells by targeting IGF-1 and regulate breast and liver cancer processes by inhibiting the expression of Cyclin B1 and IGF-1R (Khan et al. 2013;Li et al. 2017a). Studies have shown that miR-379-5p expression was significantly downregulated in cancer tissues, such as liver, bladder, and breast cancer tissues, and closely related to tumor stage; moreover, miR-379-5p can also inhibit the transformation and metastasis of endothelium to mesenchymal cells in liver cancer (Chen et al. 2016;Wu et al. 2017a). To date, miR-379 has been widely studied in human cancer, and its research on reproductive organ development has not been reported yet. WNT8A is a member of the WNT ligand family, which participates in the WNT, Hedgehog, and Hippo signaling pathway. Most studies reported that WNT8A gene plays a vital role on embryonic development and early organ formation (Narayanan 2012;Vendrell et al. 2013;Wylie et al. 2014;Hino et al. 2018). It is reported that the WNT pathway is involved in the major events during testis development, including primordial germ cell specification, proliferation and migration, testis determination, spermatogenesis, and somatic cell regulation (Dong et al. 2015). It was also reported that the WNT signaling pathway is involved in male reproductive organ development (Sun and Wang 2003), the Hedgehog signaling pathway is involved in testicular cord and spermatogonial stem cell development (Sahin et al. 2014;Li et al. 2017b), the Hippo signaling pathway is involved in regulating growth, proliferation, and apoptosis of cells, and regulating proliferation and maintenance of stem cells to control organ size (Pan 2007;Zhang et al. 2008). Therefore, it can be inferred that oar-miR-379-5p may be involved in Small Tail Han Sheep testis development and spermatogenesis by downregulated WNT8A, which involved in these pathways.
Currently, the dual luciferase reporter gene method is the most convenient method for verifying the relationship between miRNAs and their target mRNAs. The results demonstrated that oar-miR-379-5p can decrease WNT8A expression by targeting WNT8A-39UTR. It is the first time to verify that WNT8A is the target gene of oar-miR-379-5p. However, the hypothesis that oar-miR-379-5p is involved in testicular development and spermatogenesis by targeting WNT8A needs to be further studied and verified. The results of this study will provide theory evidence for identifying the molecular markers in Small Tail Han Sheep testis development and further finding out the molecular markers affecting male sheep reproductive ability.

CONCLUSION
This study is the first integrative analysis of miRNA and mRNA expression profiles in Small Tail Han Sheep testis development. We identified the key miRNAs and miRNA-targets in 2vs. 6-, 6-vs 12-, and 2-vs 12-month-old testes, which will provide theoretical basis for elucidating the molecular regulatory mechanisms underlying sheep testis development and identifying the biomarkers that affect reproductive efficiency of male Small Tail Han Sheep. Our study is the first to reveal that WNT8A was targeted by oar-miR-379-5p. We present a hypothesis that oar-miR-379-5p is involved in sheep testis development and spermatogenesis by targeting WNT8A and we will further investigate and verify this hypothesis.

ACKNOWLEDGMENTS
This work was supported by Jilin Province Science and Technology Development Project (20160204018NY), Jilin Province modern sheep industry technology system project, Jilin Province Department of Education"13th Five-Year" Science and Technology Development Project (JJKH20180691). We appreciate the Jilin Sheng Hua Animal Husbandry Co. Ltd for providing the experimental animals and Shanghai Novel Bio-Pharm Technology Co. Ltd for help with the RNA-sequencing.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE
Experimental animals were provided Songyuan Sheng Hua Animal Husbandry Co. Ltd (Jilin, China). The methods were carried out in accordance with relevant guidelines from the Ministry of Agriculture of the People's Republic of China. All experimental protocols were approved by the Jilin Laboratory Animal Specialized Committee.
AUTHORS' CONTRIBUTIONS MB and HZJ conceived and designed the experiments. MB, LMS, CJ, JRL, YH, and HL performed the experiments. MB analyzed the data, and wrote the paper. YC gave some pertinent suggestions on the experiment and writing. All authors read and approved the final manuscript.

COMPETING INTERESTS
The authors declare that there are no competing interests.