Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries) provides novel insights into the genetic mechanisms for female reproductive success

Heredity https://doi.org/10.1038/s41437-018-0090-1 ARTICLE Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries) provides novel insights into the genetic mechanisms for female reproductive success Ji Yang1●Xin Li1,2 ●Yin-Hong Cao1,2 ●Kisun Pokharel 3●Xiao-Ju Hu1,2 ●Ze-Hui Chen1,2 ●Song-Song Xu1,2 ● Jaana Peippo3●Mervi Honkatukia3●Juha Kantanen3●Meng-Hua Li1 Received: 12 January 2018 / Revised: 20 March 2018 / Accepted: 18 April 2018 © The Author(s) 2018. This article is published with open access Abstract Prolific breeds of domestic sheep (Ovis aries) are important genetic resources due to their reproductive performance, which is characterized by multiple lambs per birth and out-of-season breeding. However, the lack of a comprehensive understanding of the genetic mechanisms underlying the important reproductive traits, particularly from the evolutionary genomics perspective, has impeded the efficient advancement of sheep breeding. Here, for the first time, by performing RNA-sequencing we built a de novo transcriptome assembly of ovarian and endometrial tissues in European mouflon (Ovis musimon) and performed an mRNA–miRNA integrated expression profiling analysis of the wild species and a highly prolific domestic sheep breed, the Finnsheep. We identified several novel genes with differentially expressed mRNAs (e.g., EREG, INHBA,SPP1,AMH,TDRD5, and ZP2) between the wild and domestic sheep, which are functionally involved in oocyte and follicle development and fertilization, and are significantly (adjusted P-value < 0.05) enriched in the Gene Ontology (GO) terms of various reproductive process, including the regulation of fertilization, oogenesis, ovarian follicle development, and sperm–egg recognition. Additionally, we characterized 58 differentially expressed miRNAs and 210 associated target genes that are essential for the regulation of female reproduction cycles through specific regulatory networks [e.g., (miR-136, miR-374a,miR-9-5p)-(EREG,INHBA)]. Furthermore, our integrated mRNA and miRNA expression profiling analysis elucidated novel direct and indirect miRNA/mRNA causal regulatory relationships related to the reproductive traits of the Ovis species. This study provides in-depth insights into the genomic evolution underlying the reproductive traits of the Ovis species and valuable resources for ovine genomics. Introduction Reproduction is an important but complex biological process in mammals. From the physiological perspective, the reproductive performance of a female is attributable to a series of biological events occurring between oocyte maturation and birth (Vinet et al. 2012). An increased number of offspring is largely associated with increased ovulation rates during the estrous cycle and subsequent successful embryonic and fetal development (Vinet et al. 2012; Warriach et al. 2015). The ovary and uterus are two essential parts of the female reproductive tract. Function of the ovary includes the production of oocytes, excretion of reproductive hormones, and maintenance of female estrus cycles (Peters and McNatty 1980). The endometrium, which is the inner-most epithelial layer of the mammalian uterus, functions to prevent adhesions between the opposing walls of the myometrium, thereby maintaining the patency These authors contributed equally: Ji Yang, Xin Li, Yin-Hong Cao, Kisun Pokharel, Xiao-Ju Hu, Ze-Hui Chen. *Meng-Hua Li [email protected] 1CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences (CAS), Beijing 100101, China 2College of Life Sciences, University of Chinese Academy of Sciences (UCAS), Beijing 100049, China 3Production Systems, Natural Resources Institute Finland (Luke), 31600 Jokioinen, Finland Electronic supplementary material The online version of this article (https://doi.org/10.1038/s41437-018-0090-1) contains supplementary material, which is available to authorized users. 1234567890();,: 1234567890();,: of the uterine cavity (Fazleabas and Strakova 2002). Therefore, the ovary, endometrium, and their communications can directly affect female reproductive success (Peters and McNatty 1980; Fazleabas and Strakova 2002; Wira et al. 2015). From the evolutionary perspective, the reproductive success of females is optimally balanced by the frequency of the estrous cycle per year, time between parturition and first ovulation, number and size of offspring, and time at which sexual maturity is reached (Walker et al. 2008; Bronson 2009; Barrett et al. 2009). During and after domestication, livestock reproductive traits have been continuously selected. Wild sheep such as Asian mouflon Ovis orientalis and European mouflon Ovis musimon typically show pronounced seasonal estrous in all reproductive parameters and are typically monotonous (Lincoln 1989; Santiago-Morenoa et al. 2000; Garel et al. 2005). After domestication, many sheep (Ovis aries) breeds continue to maintain seasonal estrus and have one lamb per litter, while several prolific breeds typically demonstrate the capability of non-seasonal breeding and have multiple lambs per litter, which could be ascribed to long-term artificial and natural selection (Rasali et al. 2006). Additionally, prolific breeds of domestic sheep show much earlier sexual maturation (as early as 6 months age) than wild sheep (more than 1.5 years age; Cugnasse et al. 1985; Garel et al. 2005). Early genetic mapping studies have detected a few quantitative trait loci (QTLs) and functional genes (e.g., BMPR1B,BMP15, and GDF9) associated with high prolificacy in various sheep breeds (Hanrahan et al. 2004; Bodin et al. 2007; Abdoli et al. 2013). However, to our knowledge no study has explored the genetic mechanisms underlying the development of reproductive traits from wild to domestic sheep. Transcriptome profiling using next-generation sequencing has been successfully used to measure the levels of transcripts and their isoforms (Wang et al. 2009). In livestock, mRNA profiling has been used to characterize largescale gene expression patterns and identify differentially expressed genes that are specifically associated with production traits (Mortazavi et al. 2008; Jiang et al. 2014; Suravajhala et al. 2016; Clark et al. 2017). Additionally, microRNA (miRNA) profiling has successfully used to identify particular miRNAs, predict their target genes, and analyze their potential biological functions (Ambros 2004; Bartel 2009; McDaneld 2009; Liu et al. 2014). In sheep, previous RNA profiling investigations have revealed novel mRNAs, miRNAs, relevant functional genes, and variants associated with a variety of traits, including high prolificacy traits such as multiple births and non-seasonal breeding (Di et al. 2014; Pokharel et al. 2015; Hu et al. 2016). However, integrated analyses of large-scale mRNA and miRNA expression profiling in two (e.g., ovary and endometrium) or more female reproductive tissues from sheep are scarce (but see Miao and Qin 2015; Pokharel et al. 2018). Additionally, to date no full-length transcriptome has been available for any wild sheep, which has impeded studies investigating the genomic evolution and functional diversification of specific miRNA genes (Liu et al. 2008; Berezikov 2011; Ameres and Zamore 2013) between wild and domestic sheep. In this study, we performed genome-wide profiling of the mRNA and miRNA in ovarian and endometrial tissues in European mouflon and a prolific breed of domestic sheep, the Finnsheep (Fig. S1 and Table S1). In particular, for the first time we reconstructed a de novo transcriptome of wild sheep (Fig. S2 and Table S1) using the Trinity platform, which provides a novel algorithm for assembling full-length transcripts (Grabherr et al. 2011; Haas et al. 2013). The European mouflon is considered one of the close wild relatives to domestic sheep (Hiendleder et al. 2002; Lv et al. 2015). The Finnsheep, which is a native sheep breed from Finland, is well-known worldwide due to its excellent prolific performance, including its early maturity (at c. 6 months of age), multiple lamb births (typically 3–5 lambs per birth, and the recorded highest number of lambs is 9), and non-seasonal breeding (Maijala 1996; Li et al. 2011; Mullen and Hanrahan 2014). Nevertheless, early mRNA and miRNA studies investigating the prolificacy traits in Finnsheep have included very limited sample sizes, and the analyses of mRNAs and miRNAs were performed separately (Pokharel et al. 2015; Hu et al. 2016). By performing integrated and comparative analyses of the data within and between species, our aims are threefold: (i) to reveal the genetic mechanisms underlying the evolution of the reproduction traits in species of the genus Ovis; (ii) to decipher the genetic architecture of the high prolificacy traits in a domestic sheep breed; and (iii) to create novel mRNA and miRNA genomic resources associated with reproductive traits in wild and domestic sheep species. Our results provide a comprehensive understanding of the genetic regulatory network responsible for the high prolificacy traits and molecular mechanisms underlying their evolution in Ovis species, which could also be applied to other mammals. Materials and methods Sample collection, RNA isolation, and sequencing All animal handling was carried out under the license approved by the Southern Finland Animal Experiment Committee (approval No. ESAVI/5027/04.10.03/2012). From the Finnsheep, six ovarian samples (F1-O, F2-O, F3O, F4-O, F5-O, and F6-O) were collected from six ewes (F1, F2, F3, F4, F5, and F6), and two endometrial samples Ji Yang et al. (F2-E and F3-E) were taken from F2 and F3 individuals. The ewes in the experiment were mature, with an approximate age and weight of 4.25 years and 71.2 kg, respectively. The ovarian samples were collected at the follicular growth phase based on hormonal profiles, and endometrial samples were collected during early pregnant stage after ovulation (Pokharel et al. 2018). Three European mouflon ewes (M1, M2, and M3) were killed by hunters in autumn of 2013 (M1 and M2) and 2014 (M3). To obtain comparable sample sizes between the wild and domestic sheep, we collected six ovarian samples (M1-OA, M1-OB, M2-OA, M2-OB, M3-OA, and M3-OB) representing two replicates (coded as -OA and -OB) from each of the three European mouflons. Two endometrial samples (M2-EA and M2-EB) were taken by cytobrush from both uterine horns of the older individual M2 (5–6 years old), and we did not take endometrial samples from the other two younger wild sheep. We did not have any information about the stage of the cycle in the mouflons except that the older one (M2) had a corpus luteum. In total, 16 samples consisting of 12 ovarian samples and 4 endometrial samples were included in the mRNA profiling analysis. Additionally, all 16 samples were included in the miRNA profiling analysis, except for M2-EA (Table S1). After the sample collection, the tissue samples of domestic sheep were quickly stored in RNAlater reagent (Ambion/Qiagen, Valencia, CA, USA) per the manufacturer’s instructions and transported to the laboratory. For the European mouflons, the collected samples were stored in RNAlater within 1 h of death due to the forest areas. The tissue samples were then transferred to be stored at −80°C in laboratory until extraction. We used an RNeasy plus mini kit (Qiagen, Valencia, CA, USA) to extract the mRNA and miRNA from the tissues according to the manufacturer’s protocol. The RNA concentration and RNA integrity number were measured using a Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany). Libraries of mRNAs and miRNAs were prepared using the Illumina’s TruSeq library preparation kits and sequenced using Illumina Hiseq2000 at the Finnish Microarray and Sequencing Center, Turku, Finland. The mRNA libraries were sequenced using 100 base-pair (bp) paired-end sequencing chemistry, whereas the miRNAs were sequenced using the single-end 50 bp approach. Quality control and mapping FastQC v.0.11.4 (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/) was used to examine the quality of the raw RNA-Seq (both the mRNA-Seq and miRNA-Seq) reads, and the adapter sequences were removed using the program Cutadapt v.1.9.1 (Martin 2011). The obtained trimmed mRNA-Seq reads were then mapped against the sheep reference genome (Oar v.4.0) supported by the gene annotation file of Ensembl release 83 using the program TopHat v.2.0.8b (Trapnell et al. 2012). For the miRNA-Seq data, the program PRINSEQ-LITE v.0.20.4 (Schmieder and Edwards 2011) was used for further miRNA read preprocessing (including filtering for reads containing ambiguous reads, data reformatting, and adapter trimming). Clean reads between 18 and 26 nucleotides (nt) in length were used for the subsequent analysis. Then, the Bowtie-build tool was used to build an index for the genome alignment using the sheep genome assembly Oar v.3.1 (ftp://ftp. ensembl.org/pub/release-83/fasta/ovis_aries/dna/), and the clean reads (18–26 nt) were mapped to the index using the program Bowtie v.2.2.1 (Langmead et al. 2009). The bioinformatics pipeline used to analyze the mRNA-Seq and miRNA-Seq data is summarized in Fig. S1. De novo transcriptome assembly and annotation of European mouflon We used the Trinity v.2.1.1 package (Grabherr et al. 2011) to construct a de novo transcriptome assembly of European mouflon using the ovarian (M1-OA, M1-OB, M2-OA, M2OB, M3-OA, and M3-OB) and endometrial (M2-EA and M2-EB) mRNA sequences. The quality of the de novo assembly was assessed by examining the proportion of reads mapped to the assembly. We used the RSEM (RNASeq by expectation-maximization) tool (Li and Dewey 2011) to estimate the number of RNA-Seq fragments that map to each contig in the individual M2 (M2-OA, M2-OB, M2-EA, and M2-EB), and detected the differentially expressed genes between the two types of tissues, i.e., ovary and endometrium, using the EdgeR package (Robinson et al. 2010). Unigenes that had |log2(fold change)| ≥2 and a false discovery rate ≤0.01 were considered significantly differentially expressed. The differentially expressed unigenes were mapped to the NCBI nt collection database using the BlastN (Altschul et al. 1997) tool in the Blast v.2.4.0+package (http://blast.ncbi.nlm.nih.gov) with a matched standard “e-value =1e-10 and matching ratio = 100%”. The program PANTHER v.10 (Mi et al. 2016) was used to perform the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, and the threshold was set to a P-value < 0.05 after the Bonferroni correction. The bioinformatics analysis workflow is summarized in Fig. S2. mRNA differential expression and functional analysis After the mapping, the programs Cufflinks, Cuffmerge, and Cuffdiff in the Cufflinks software package v.2.1.1 (Trapnell et al. 2012) were used to assemble the mRNA clean reads Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries). . . and calculate the gene expression levels based on the fragments per kilobase of exon per million fragments mapped (FPKM) values. First, we performed the mRNA transcript assembly and created a “merged transcripts.gtf” for all ovarian or endometrial samples from one species using the programs Cufflinks and Cuffmerge. Then, the program Cuffdiff was used to calculate the level of gene expression and perform differential expression analyses between the two species (i.e., ovary/endometrium: European mouflon vs. Finnsheep) or two types of tissues (i.e., European mouflon: ovary vs. endometrium). The differentially expressed genes (DEGs) were defined as genes that showed |log2(fold change)| ≥2 and Padj (i.e. adjusted Pvalue) ≤0.01 in the comparisons. Furthermore, we performed GO and KEGG pathway enrichment analyses of functional genes associated with the differentially expressed mRNAs using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tools v.6.8 (Huang et al. 2009). The enriched GO terms were considered significant at an adjusted P-value < 0.05. A word cloud was created for the GO terms and pathways using an online Wordcloud generator (https://www.jasondavies.com/wordcloud/). miRNA quantification and differential expression analysis The miRDeep2 v.0.0.7 software package (Friedländer et al. 2012) was used to predict the putative precursor miRNA sequences, and precursor miRNAs with 10 or fewer reads and miRDeep2 scores < 5 were removed. Then, the filtered miRNAs were aligned to the miRNA database miRBase v.21 (http://www.mirbase.org/) to identify mature miRNAs that contain sheep miRNAs (i.e., miRNAs deposited in the miRBase of O. aries), conserved miRNAs (i.e., miRNAs matching those in other mammals in miRBase with fewer than four mismatches), and novel miRNAs (remaining miRNAs). The differentially expressed ovarian and endometrial miRNAs (including the sheep and conserved miRNAs) between European mouflon and Finnsheep and between the two types of tissues (ovary vs. endometrium) in European mouflon were identified using the DESeq v.1.22.0 program (Anders and Huber 2010). By applying the threshold values of Padj (i.e., adjusted P-value) ≤0.05 and |log2(fold change)| ≥1, the remaining miRNAs considered to be upregulated or downregulated in one of the two species or types of tissues were used in the following mRNA–miRNA integrated analyses. miRNA–mRNA regulatory network analysis The program TargetScan v.7.0 (Agarwal et al. 2015) was used to predict the binding of the differentially expressed miRNAs to their putative targets, i.e., the 5′seed regions (2–8 nt) of the miRNAs complementary to the conserved 8and 7-mer sites in the 3′-untranslated region of the mRNA (Lewis et al. 2005). miRNAs typically function in posttranscriptional regulation by suppressing or silencing specific target genes. Thus, we searched for inverse correlations between the expression levels of the miRNAs and their predicted target genes, i.e., between the upregulated miRNAs and downregulated mRNAs, and vice versa. To narrow the list of potential miRNA–target interactions, the putative target genes were overlapped with the DEGs in the mRNA profiling data. Then, target genes with a context++ score percentile lower than 50 were removed. Furthermore, the upregulated and downregulated miRNAs and their target genes were included in the network visualization analysis using the program Cytoscape v.3.4.0 (Shannon et al. 2003). Only the top 10 downregulated miRNAs were considered core miRNAs and used to build the main regulatory network. Furthermore, GO and pathway enrichment analyses of the miRNA target genes were conducted using DAVID tools v.6.8 (Huang et al. 2009). GO terms and pathways with an adjusted P-value < 0.05 were considered significantly enriched. A flowchart of the mRNA–miRNA integrated analysis is provided in Fig. S1. To confirm the credibility of the inferred miRNA and target gene pairs, we searched the miRTarBase (http://mirta rbase.mbc.nctu.edu.tw/; Chou et al. 2016) database and previous studies for differentially expressed miRNAs, predicted target genes, and miRNA–gene pairs consistent with those identified in this study. Results and Discussion mRNA and miRNA profiling In total, we generated approximately 127.1 Gb mRNA and 21.3 Gb miRNA raw data. After filtering, we obtained a large number of high-quality reads: 1083.6 million mRNA and 56.8 million miRNA reads in European mouflon, and 934.6 million mRNA and 127.8 million miRNA reads in Finnsheep. Most high-quality mRNA and miRNAs have been mapped to the O.aries reference genome Oar_v.4.0 (Table S2). Regarding the mRNA expression, more than 89% of the genes were expressed at <50 FPKM and only <1.1% of the genes were expressed at >500 FPKM in all four tissue-species combinations (Table S3). Altogether, 23 557 common genes were expressed in the ovaries and 19 208 common genes were expressed in the endometria (Table S3). In the de novo assembly of the European mouflon transcriptome, 1083.6 million clean mRNA reads were assembled into 844 129 unigenes, with an N50 length of 700 bp Ji Yang et al. and an average length of 562 bp (Table 1). The observed length of N50 is comparable to that reported in domestic sheep (e.g., 508–1482 bp) (Yue et al. 2015; Zhang et al. 2015), thus highlighting the high quality of the assembled European mouflon transcriptome. Based on the transcriptome, we identified 5739, 5831, 5897, and 5981 annotated DEGs in the four comparisons between the ovary and endometrium tissues (M2-OA vs. M2-EA, M2-OA vs. M2-EB, M2-OB vs. M2-EA, and M2-OB vs. M2-EB) in European mouflon, respectively. However, when the analyses were performed using the sheep reference genome (Oar v.4.0), we found far fewer DEGs (e.g., 1925 DEGs for M2-O-B vs. M2-E-B). Despite the possible implication of using the O. aries genome as a reference, the observation suggests that the de novo approach is accurate and that the quality of the de novo assembly is good. We identified 373 mature miRNAs in the ovaries and 192 mature miRNAs in the endometria, including 89 and 59 sheep miRNAs, 174 and 104 conserved miRNAs, and 110 and 29 novel miRNAs in the two tissues, respectively (Table S4). High levels of expression were observed in the sheep miRNAs, followed by conserved miRNAs and novel miRNAs (Table S4). Of the 122 novel miRNAs (Table S5), 20 miRNAs were also identified as novel miRNAs and validated in our previous study (Hu et al. 2016). The chromosomal distribution of the sheep and conserved miRNAs indicated the presence of large numbers of miRNAs on chromosomes 18 and X, whereas none of the miRNAs were expressed on chromosome 8 based on the ovarian data and chromosomes 8, 17, and 25 based on the endometrium data (Table S6). mRNAs and miRNAs showing the greatest expression Of the genes exhibiting the greatest expression (FPKM > 3000) in ovary, 7 genes in the Finnsheep and 15 genes in the European mouflon encode ribosomal protein (Fig. 1a and Table S7). In the endometrium, only 1 top gene identified in the Finnsheep and 22 top genes identified in the European mouflon encode ribosomal protein involved in the ribosome pathway (Fig. 1b and Table S7). The ribosome pathway has been reported to be among the most significant pathways involved in early embryogenesis (Berres et al. 2017). The observed large number of ribosome family genes in both tissues of European mouflon could be one explanation for the low reproductive success of the wild sheep (i.e., seasonal estrus, 1.5 years of age to sexual maturation and one lamb per litter; Cugnasse et al. 1985; Garel et al. 2005) compared to the prolificacy of domestic sheep, including Finnsheep (i.e., out-of-season breeding, 6 months of age to sexual maturation and three lambs per litter on average; Maijala 1996; Li et al. 2011; Mullen and Hanrahan 2014). Moreover, the observed smaller number of ribosome family genes expressed at a higher level in Finnsheep than in European mouflon is likely due to the evolution of reproductive traits between the two species during the domestication and long-term artificial selection processes. Additionally, this observation could be related to the environment of the mouflon being harsher than the domestic environment of the Finnsheep; in such a harsh environment, it may be beneficial to have fewer offspring. Of the 10 most highly expressed ovarian miRNAs, 7 (miR-10b,miR-143,miR-26a,miR-27b,let-7f,miR-191, and miR-22-3p) miRNAs are common between European mouflon and Finnsheep (Table S8). Of the endometrial miRNAs, three (miR-10b,miR-100, and miR-34c-5p) miRNAs are highly expressed in both species (Table S8). We also observed tissueand species-specific miRNAs (Table S9), which could possibly be associated with particular functions. For instance, miR-181a, which is specifically and abundantly expressed in the Finnsheep ovary, was demonstrated to be associated with pregnancy complications and was more highly expressed in normal term than in pre-term delivery (Mayor-Lynn et al. 2011). Additionally, miR-186, which is specifically and abundantly expressed in the European mouflon endometrium, was observed to be more highly expressed in endometrial cancer (Zhou et al. 2008). Differentially expressed mRNAs and miRNAs In total, 192 and 1065 DEGs responsible for ovarian and endometrial mRNAs were identified between the two species, respectively. Of these DEGs, 149 ovarian and 492 endometrial DEGs were upregulated in European mouflon, Table 1 Summary of constituent data of Trinity-assembled European mouflon transcriptome Description Number Before trimming Raw reads 1 083 877 020 After trimming Clean reads 1 083 600 777 Average length of clean reads (bp) 101 N percentage 0 GC content (%) 48.92 After assembly Total Trinity “genes”844 129 Total Trinity transcripts 1 102 841 Statistics based on all transcript contigs Mean length of contigs (bp) 947.39 N50 (bp) of contigs 2169 Statistics based on unigenes Mean length of unigenes (bp) 562.14 N50 (bp) of unigenes 700 Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries). . . while 43 ovarian and 573 endometrial DEGs were upregulated in Finnsheep (Tables S10 and S11). Additionally, we identified 53 and 5 differentially expressed ovarian and endometrial miRNAs between the two species, respectively (Table 2). Based on the differentially expressed mRNAs and miRNAs in the tissues, heat maps were constructed (Fig. 2) and showed distinct expression patterns between European mouflon and Finnsheep, indicating that the evolution of the geneand post-transcriptional regulation could have contributed to the improved reproductive performance during the domestication process. Of the upregulated genes in the ovary and endometrium tissues from the Finnsheep, 4 genes (CYP19,PTGS2,FST, and PAPSS2) are associated with early maturity, 14 genes (EREG,INHBA,SERPINA5,SPP1,LIPG,PTGS2,AREG, TDGF1,FST,C4BPA,SERPINE2,TRPV3,LDLR, and TINAGL1) are associated with fecundity, and two genes (CYP19 and SERPINE2) are associated with non-seasonal breeding. In addition, we compared the 192 ovarian DEGs with the QTLs for the reproductive traits in the Sheep QTL database (https://www.animalgenome.org/cgi-bin/QTLdb/ OA/), and found 8 genes (NME9,LRRIQ4,SPATC1L, PRDM15,PEX5L,COL6A6,RARRES1, and SPOCK1; Table S12) overlapped with previously identified QTLs associated with traits such as reproductive seasonality and total lambs born. Similarly, we found 73 endometrial DEGs that overlapped with the reported QTLs for the reproductive seasonality, and 1 gene, i.e., RASL11B, was overlapped with the QTL chromosome 6 as follows: 68295117–68445447 for the total lambs born trait (Table S13). GO and KEGG enrichments of differentially expressed mRNAs In the functional enrichment analysis of the 192 ovarian DEGs between European mouflon and Finnsheep, the top enriched GO categories were mainly relevant to reproduction, including fertilization, reproductive developmental process, reproductive process, oogenesis, and ovarian follicle development (Fig. S3a and Table S14). Fertilization refers to the fusion of gametes that initiates the development of a new individual organism (Merriam-Webster 2017). During the process, the ZP proteins (ZP1 to ZP4) in the zona pellucida function to surround the egg, bind the sperm and mediate species-selective sperm–oocyte interactions (Suzuki et al. 2015; Avella et al. 2016). ZP2 peptide beads have been shown to markedly inhibit the fertilization of ovulated eggs, prevent sperm binding in mice and humans (Gupta et al. 2012; Avella et al. 2016) and induce infertility in female mice (Avella et al. 2016). ASTL, which is involved in the negative regulation of fertilization by the ZP2 and ZP4 genes, plays a definitive role in ensuring monospermic fertilization by encoding ovastacin, which is a cortical granule protease (Burkart et al. 2012). Thus, the significantly upregulated expression of the ASTL,ZP2,ZP3, and ZP4 genes in the European mouflon ovary (Table S10) may indicate that these genes play a role in the low fertility of wild sheep compared to the high prolificacy of domestic sheep (Fig. 3). The GO terms of the reproductive development process and ovarian follicle development included the genes EREG, INHBA,SERPINA5, and SPP1, which are significantly upregulated in Finnsheep (Tables S10 and S14). EREG is a member of the epidermal growth factor family, which plays an important role in mediating the luteinizing hormoneinduced ovulation process in preovulatory follicles (Amsterdam 2010; Kim et al. 2011). INHBA plays an essential role in early oocyte and follicular development (Thomas et al. 2003) and has a significant effect on the litter size in sheep (Leyhe et al. 1994) and goat (Zi et al. 2012). SERPINA5, which encodes a serine protease inhibitor, may be involved in follicular growth and atresia in the bovine Fig. 1 Number of genes exhibiting the greatest expression (FPKM > 3000) in the Finnsheep and European mouflon ovaries (a) and endometria (b), and their overlap with genes in the ribosome pathway Ji Yang et al. ovarian follicle (Hayashi et al. 2011). SPP1 is located on sheep chromosome 6 close to the high prolificacy gene FecB in sheep (Montgomery et al. 1993; Montgomery et al. 1994). Additionally, the genes AMH,TDRD5,SOHLH1, and NOBOX are upregulated in European mouflon (Tables S10 and S14). AMH is a central influencing factor in gonad development, and could inhibit germ cell proliferation and differentiation in both sexes (Durlinger et al. 2002). Recent clinical data have indicated that serum AMH level is a useful indicator of ovarian aging and premature ovarian failure in women (de Vet et al. 2002; Alipour et al. 2015). The TDRD5 protein is an evolutionarily conserved protein involved in germ cell development. Injection of TDRD5deficient round spermatids into oocytes result in fertile offspring (Yabuta et al. 2011). SOHLH1 is a critical regulator of oogenesis and can downregulate NOBOX, which functions to disrupt folliculogenesis (Pangas et al. 2006). Therefore, our results suggest that the upregulated Table 2 List of significantly differentially expressed miRNAs between Finnsheep and European mouflon tissues (ovary and endometrium) miRNA log2(fold change)aPadjbmiRNA log2(fold change)aPadjb Ovary (European mouflon up or Finnsheep down) Ovary (European mouflon down or Finnsheep up) miR-136 7.70 2.03E-40 miR-140-3p 4.22 2.60E-30 miR-708-3p 6.01 1.40E-08 miR-370-5p 3.95 1.88E-13 miR-101 4.62 4.22E-05 miR-615-3p 3.87 3.09E-08 miR-374a 4.54 1.52E-15 miR-197 3.50 3.28E-03 miR-148a 4.47 1.07E-14 miR-1306 3.36 9.62E-14 miR-3959-5p 4.06 1.19E-16 miR-1247-5p 3.20 5.09E-03 miR-411a-5p 3.56 9.57E-27 miR-744 3.14 2.94E-11 miR-6119-5p 3.10 6.12E-03 miR-296-3p 2.94 1.75E-17 miR-335 3.03 4.52E-03 miR-432 2.54 1.10E-04 miR-9-5p 2.73 8.40E-03 miR-30d 2.54 1.92E-22 miR-369-5p 2.45 3.26E-05 miR-423-3p 2.39 3.17E-13 miR-340-5p 2.41 1.24E-09 miR-331 2.28 8.14E-03 miR-21 2.24 9.93E-20 miR-2284x 2.24 1.59E-13 miR-27a 2.13 3.16E-03 miR-193b-3p 1.98 6.84E-03 miR-30e-5p 2.07 9.30E-11 miR-16b 1.81 1.97E-12 miR-143 1.99 8.11E-03 miR-665-5p 1.71 1.60E-04 miR-99a 1.96 8.07E-26 miR-485-5p 1.67 1.62E-03 miR-374b 1.90 2.14E-03 miR-30f 1.65 5.96E-05 miR-10b 1.81 8.11E-03 let-7e-5p 1.64 2.69E-07 miR-30b 1.71 3.39E-02 miR-382-5p 1.47 9.02E-05 miR-379-5p 1.38 1.64E-09 miR-92a 1.39 5.07E-08 miR-98 1.38 1.45E-02 miR-28-3p 1.39 2.15E-05 miR-27b 1.20 7.73E-05 miR-378 1.31 7.16E-03 miR-409-5p 1.28 2.55E-03 miR-15b 1.26 1.27E-05 miR-3432-5p 1.20 3.95E-02 let-7c 1.19 1.40E-08 miR-93 1.18 3.06E-04 miR-22-3p 1.17 5.07E-08 miR-874 1.08 2.00E-02 Endometrium (European mouflon up or Finnsheep down) Endometrium (European mouflon down or Finnsheep up) miR-34c-5p 4.73 9.99E-03 miR-378 3.18 2.84E-02 miR-10b 3.82 4.92E-07 miR-100 3.29 2.12E-03 miR-92a 2.18 2.39E-02 aThe threshold for fold change is set to log2(fold change) > 1 or <−1 bThe statistical significance is evaluated as Padj (i.e., adjusted P-value) < 0.05 Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries). . . expression of these genes (i.e., EREG,INHBA,SERPINA5, and SPP1 in Finnsheep, and AMH,SOHLH1,TDRD5, and NOBOX in European mouflon) is most likely associated with the hyper-prolificacy of domestic sheep (e.g., Finnsheep) and the single lambing of wild sheep (e.g., European mouflon; Fig. 3). We also observed other top GO categories Fig. 2 Heat maps of differentially expressed mRNAs and miRNAs between Finnsheep and European mouflon. aHeat map showing the expression intensity of 1065 differentially expressed genes (DEGs) based on a comparison of endometria from the two species. bHeat map showing the expression intensity of 192 DEGs based on a comparison of ovaries from the two species. c Heat map showing the expression intensity of 5 differentially expressed miRNAs based on a comparison of endometria from the two species. dHeat map showing the expression intensity of 53 differentially expressed miRNAs based on a comparison of ovaries from the two species. Heat maps were generated using the R heat map package v.1.0.8 (http://cran.r-project.org/web/pa ckages/pheatmap/), which is based on the k-means clustering algorithm and the Boolean values to determine if the rows (i.e., mRNAs or miRNAs) should be clustered. The important genes associated with the reproductive traits such as early mature, fecundity, nonseasonal breeding, and total lambs born are labeled in the figure Ji Yang et al. related to immune processes such as the humoral immune response and humoral immune response mediated by circulating immunoglobulin (Table S14). The KEGG analysis of the 192 ovarian DEGs identified three significant (adjusted P-value < 0.05) pathways, i.e., ovarian steroidogenesis, complement and coagulation cascades, and MAPK signaling pathways (Table S14). The MAPK (i.e., p38αMAPK) pathway has been shown to play an essential role in female reproduction, particularly in follicular development (Hu et al. 2017) (Fig. 3). In the functional enrichment of the endometrial DEGs between the wild and domestic sheep, the GO terms of the 492 upregulated genes (RARRES2,LTBP1, and TGFB3)in European mouflon were most significantly enriched in the developmental process, such as anatomical structure development (Fig. S3b and Table 3), whereas the top enriched GO terms of the 573 upregulated genes (PTGS2, LDLR, and SREBF1) in Finnsheep included biosynthetic and metabolic processes, such as lipid metabolic process (Fig. S3b and Table 3). An early study investigating the bovine endometrial transcriptome indicated an earlier shift from proliferation to metabolism at early diestrus in the cow (Mesquita et al. 2015). Cell proliferation was enriched in the small follicle endometrium, whereas biosynthetic and metabolic processes were enriched in the large follicle endometrium, which presented higher proliferative activity in the luminal epithelium, glandular epithelium, and stroma than in the small follicle endometrium (Mesquita et al. 2015). Additionally, the large follicle endometrium is an active endometrial phenotype that is associated with the optimal uterine environment (Miller and Moore 1976; Mesquita et al. 2015). Thus, Finnsheep could have a larger follicle endometrium and better receptivity of the uterus than European mouflon (Fig. 3). In the KEGG analysis of the 1065 endometrial DEGs, the top enriched pathways were enriched in the extracellular matrix (ECM)–receptor interaction, steroid biosynthesis, and biosynthetic and metabolic processes (Table S15). The ECM, which is a complex matrix of biological macromolecules, can transfer signals to cells via surface receptors during cell adhesion and other various signaling transduction pathways. Additionally, the ECM can send signals to the cytoplasm and nucleus to influence gene expression or cellular activities (Alldinger et al. 2006). Thus, adhesion and connection mechanisms could play important roles in the endometrium at the molecular level, and different cell adhesion may occur between wild and domestic sheep. In addition, the differentially expressed mRNAs between the ovarian and endometrial tissues from the European mouflon and their functional analysis based on the reference genome and de novo approaches are presented in the Supplementary Results and Tables S16–S18. Post-transcriptional regulatory network of miRNAs and target genes We obtained 100, 64, 17, and 51 target genes for the upand downregulated miRNAs in particular tissues (Table 4). In the main downregulated network of the Finnsheep ovary, which comprised 10 downregulated miRNAs and 57 upregulated target genes (Fig. 4a), several important miRNAs Fig. 3 Functional roles of the main mRNAs and miRNAs in the endometrium and ovary of Ovis species. Ovarian and endometrial network of differentially expressed mRNAs and top 10 Finnsheep downregulated (i.e., European mouflon upregulated) miRNAs and target genes are shown, along with their functional roles associated with reproduction (bold type). Regulation and alteration of these functional roles are responsible for the evolution of high prolificacy traits. Red and blue arrows indicate downregulation and upregulation of the reproductive organs, respectively Comparative mRNA and miRNA expression in European mouflon (Ovis musimon) and sheep (Ovis aries). . .