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

He edi y h ps://doi.o g/10.1038/s41437-018-0090-1 ARTICLE Compa a i e mRNA and miRNA exp ession in Eu opean mouflon (O is musimon) and sheep (O is a ies) p o ides no el insigh s in o he gene ic mechanisms o emale ep oduc i e success Ji Yang1●Xin Li1,2 ●Yin-Hong Cao1,2 ●Kisun Pokha el 3●Xiao-Ju Hu1,2 ●Ze-Hui Chen1,2 ●Song-Song Xu1,2 ● Jaana Peippo3●Me i Honka ukia3●Juha Kan anen3●Meng-Hua Li1 Recei ed: 12 Janua y 2018 / Re ised: 20 Ma ch 2018 / Accep ed: 18 Ap il 2018 © The Au ho (s) 2018. This a icle is published wi h open access Abs ac P olific b eeds o domes ic sheep (O is a ies) a e impo an gene ic esou ces due o hei ep oduc i e pe o mance, which is cha ac e ized by mul iple lambs pe bi h and ou -o -season b eeding. Howe e , he lack o a comp ehensi e unde s anding o he gene ic mechanisms unde lying he impo an ep oduc i e ai s, pa icula ly om he e olu iona y genomics pe spec i e, has impeded he e ficien ad ancemen o sheep b eeding. He e, o he fi s ime, by pe o ming RNA-sequencing we buil a de no o ansc ip ome assembly o o a ian and endome ial issues in Eu opean mouflon (O is musimon) and pe o med an mRNA–miRNA in eg a ed exp ession p ofiling analysis o he wild species and a highly p olific domes ic sheep b eed, he Finnsheep. We iden ified se e al no el genes wi h di e en ially exp essed mRNAs (e.g., EREG, INHBA,SPP1,AMH,TDRD5, and ZP2) be ween he wild and domes ic sheep, which a e unc ionally in ol ed in oocy e and ollicle de elopmen and e iliza ion, and a e significan ly (adjus ed P- alue < 0.05) en iched in he Gene On ology (GO) e ms o a ious ep oduc i e p ocess, including he egula ion o e iliza ion, oogenesis, o a ian ollicle de elopmen , and spe m–egg ecogni ion. Addi ionally, we cha ac e ized 58 di e en ially exp essed miRNAs and 210 associa ed a ge genes ha a e essen ial o he egula ion o emale ep oduc ion cycles h ough specific egula o y ne wo ks [e.g., (miR-136, miR-374a,miR-9-5p)-(EREG,INHBA)]. Fu he mo e, ou in eg a ed mRNA and miRNA exp ession p ofiling analysis elucida ed no el di ec and indi ec miRNA/mRNA causal egula o y ela ionships ela ed o he ep oduc i e ai s o he O is species. This s udy p o ides in-dep h insigh s in o he genomic e olu ion unde lying he ep oduc i e ai s o he O is species and aluable esou ces o o ine genomics. In oduc ion Rep oduc ion is an impo an bu complex biological p o- cess in mammals. F om he physiological pe spec i e, he ep oduc i e pe o mance o a emale is a ibu able o a se ies o biological e en s occu ing be ween oocy e ma u a ion and bi h (Vine e al. 2012). An inc eased numbe o o sp ing is la gely associa ed wi h inc eased o ula ion a es du ing he es ous cycle and subsequen success ul emb yonic and e al de elopmen (Vine e al. 2012; Wa iach e al. 2015). The o a y and u e us a e wo essen ial pa s o he emale ep oduc i e ac . Func ion o he o a y includes he p oduc ion o oocy es, exc e ion o ep oduc i e ho mones, and main enance o emale es us cycles (Pe e s and McNa y 1980). The endome ium, which is he inne -mos epi helial laye o he mammalian u e us, unc ions o p e en adhesions be ween he opposing walls o he myome ium, he eby main aining he pa ency These au ho s con ibu ed equally: Ji Yang, Xin Li, Yin-Hong Cao, Kisun Pokha el, Xiao-Ju Hu, Ze-Hui Chen. *Meng-Hua Li [email p o ec ed] 1CAS Key Labo a o y o Animal Ecology and Conse a ion Biology, Ins i u e o Zoology, Chinese Academy o Sciences (CAS), Beijing 100101, China 2College o Li e Sciences, Uni e si y o Chinese Academy o Sciences (UCAS), Beijing 100049, China 3P oduc ion Sys ems, Na u al Resou ces Ins i u e Finland (Luke), 31600 Jokioinen, Finland Elec onic supplemen a y ma e ial The online e sion o his a icle (h ps://doi.o g/10.1038/s41437-018-0090-1) con ains supplemen a y ma e ial, which is a ailable o au ho ized use s. 1234567890();,: 1234567890();,: o he u e ine ca i y (Fazleabas and S ako a 2002). The e o e, he o a y, endome ium, and hei communica- ions can di ec ly a ec emale ep oduc i e success (Pe e s and McNa y 1980; Fazleabas and S ako a 2002; Wi a e al. 2015). F om he e olu iona y pe spec i e, he ep o- duc i e success o emales is op imally balanced by he equency o he es ous cycle pe yea , ime be ween pa - u i ion and fi s o ula ion, numbe and size o o sp ing, and ime a which sexual ma u i y is eached (Walke e al. 2008; B onson 2009; Ba e e al. 2009). Du ing and a e domes ica ion, li es ock ep oduc i e ai s ha e been con inuously selec ed. Wild sheep such as Asian mouflon O is o ien alis and Eu opean mouflon O is musimon ypically show p onounced seasonal es ous in all ep oduc i e pa ame e s and a e ypically mono onous (Lincoln 1989; San iago-Mo enoa e al. 2000; Ga el e al. 2005). A e domes ica ion, many sheep (O is a ies) b eeds con inue o main ain seasonal es us and ha e one lamb pe li e , while se e al p olific b eeds ypically demons a e he capabili y o non-seasonal b eeding and ha e mul iple lambs pe li e , which could be asc ibed o long- e m a ificial and na u al selec ion (Rasali e al. 2006). Addi- ionally, p olific b eeds o domes ic sheep show much ea lie sexual ma u a ion (as ea ly as 6 mon hs age) han wild sheep (mo e han 1.5 yea s age; Cugnasse e al. 1985; Ga el e al. 2005). Ea ly gene ic mapping s udies ha e de ec ed a ew quan i a i e ai loci (QTLs) and unc ional genes (e.g., BMPR1B,BMP15, and GDF9) associa ed wi h high p olificacy in a ious sheep b eeds (Han ahan e al. 2004; Bodin e al. 2007; Abdoli e al. 2013). Howe e , o ou knowledge no s udy has explo ed he gene ic mechan- isms unde lying he de elopmen o ep oduc i e ai s om wild o domes ic sheep. T ansc ip ome p ofiling using nex -gene a ion sequen- cing has been success ully used o measu e he le els o ansc ip s and hei iso o ms (Wang e al. 2009). In li e- s ock, mRNA p ofiling has been used o cha ac e ize la ge- scale gene exp ession pa e ns and iden i y di e en ially exp essed genes ha a e specifically associa ed wi h p o- duc ion ai s (Mo aza i e al. 2008; Jiang e al. 2014; Su a ajhala e al. 2016; Cla k e al. 2017). Addi ionally, mic oRNA (miRNA) p ofiling has success ully used o iden i y pa icula miRNAs, p edic hei a ge genes, and analyze hei po en ial biological unc ions (Amb os 2004; Ba el 2009; McDaneld 2009; Liu e al. 2014). In sheep, p e ious RNA p ofiling in es iga ions ha e e ealed no el mRNAs, miRNAs, ele an unc ional genes, and a ian s associa ed wi h a a ie y o ai s, including high p olificacy ai s such as mul iple bi hs and non-seasonal b eeding (Di e al. 2014; Pokha el e al. 2015; Hu e al. 2016). Howe e , in eg a ed analyses o la ge-scale mRNA and miRNA exp ession p ofiling in wo (e.g., o a y and endome ium) o mo e emale ep oduc i e issues om sheep a e sca ce (bu see Miao and Qin 2015; Pokha el e al. 2018). Addi- ionally, o da e no ull-leng h ansc ip ome has been a ailable o any wild sheep, which has impeded s udies in es iga ing he genomic e olu ion and unc ional di e - sifica ion o specific miRNA genes (Liu e al. 2008; Be e- ziko 2011; Ame es and Zamo e 2013) be ween wild and domes ic sheep. In his s udy, we pe o med genome-wide p ofiling o he mRNA and miRNA in o a ian and endome ial issues in Eu opean mouflon and a p olific b eed o domes ic sheep, he Finnsheep (Fig. S1 and Table S1). In pa icula , o he fi s ime we econs uc ed a de no o ansc ip ome o wild sheep (Fig. S2 and Table S1) using he T ini y pla o m, which p o ides a no el algo i hm o assembling ull-leng h ansc ip s (G abhe e al. 2011; Haas e al. 2013). The Eu opean mouflon is conside ed one o he close wild ela i es o domes ic sheep (Hiendlede e al. 2002; L e al. 2015). The Finnsheep, which is a na i e sheep b eed om Finland, is well-known wo ldwide due o i s excellen p olific pe o mance, including i s ea ly ma u i y (a c. 6 mon hs o age), mul iple lamb bi hs ( ypically 3–5 lambs pe bi h, and he eco ded highes numbe o lambs is 9), and non-seasonal b eeding (Maijala 1996; Li e al. 2011; Mullen and Han ahan 2014). Ne e heless, ea ly mRNA and miRNA s udies in es iga ing he p olificacy ai s in Finnsheep ha e included e y limi ed sample sizes, and he analyses o mRNAs and miRNAs we e pe o med sepa- a ely (Pokha el e al. 2015; Hu e al. 2016). By pe o ming in eg a ed and compa a i e analyses o he da a wi hin and be ween species, ou aims a e h ee old: (i) o e eal he gene ic mechanisms unde lying he e olu ion o he ep o- duc ion ai s in species o he genus O is; (ii) o deciphe he gene ic a chi ec u e o he high p olificacy ai s in a domes ic sheep b eed; and (iii) o c ea e no el mRNA and miRNA genomic esou ces associa ed wi h ep oduc i e ai s in wild and domes ic sheep species. Ou esul s p o- ide a comp ehensi e unde s anding o he gene ic eg- ula o y ne wo k esponsible o he high p olificacy ai s and molecula mechanisms unde lying hei e olu ion in O is species, which could also be applied o o he mammals. Ma e ials and me hods Sample collec ion, RNA isola ion, and sequencing All animal handling was ca ied ou unde he license app o ed by he Sou he n Finland Animal Expe imen Commi ee (app o al No. ESAVI/5027/04.10.03/2012). F om he Finnsheep, six o a ian samples (F1-O, F2-O, F3- O, F4-O, F5-O, and F6-O) we e collec ed om six ewes (F1, F2, F3, F4, F5, and F6), and wo endome ial samples Ji Yang e al. (F2-E and F3-E) we e aken om F2 and F3 indi iduals. The ewes in he expe imen we e ma u e, wi h an app ox- ima e age and weigh o 4.25 yea s and 71.2 kg, espec- i ely. The o a ian samples we e collec ed a he ollicula g ow h phase based on ho monal p ofiles, and endome ial samples we e collec ed du ing ea ly p egnan s age a e o ula ion (Pokha el e al. 2018). Th ee Eu opean mouflon ewes (M1, M2, and M3) we e killed by hun e s in au umn o 2013 (M1 and M2) and 2014 (M3). To ob ain compa - able sample sizes be ween he wild and domes ic sheep, we collec ed six o a ian samples (M1-OA, M1-OB, M2-OA, M2-OB, M3-OA, and M3-OB) ep esen ing wo eplica es (coded as -OA and -OB) om each o he h ee Eu opean mouflons. Two endome ial samples (M2-EA and M2-EB) we e aken by cy ob ush om bo h u e ine ho ns o he olde indi idual M2 (5–6 yea s old), and we did no ake endome ial samples om he o he wo younge wild sheep. We did no ha e any in o ma ion abou he s age o he cycle in he mouflons excep ha he olde one (M2) had a co pus lu eum. In o al, 16 samples consis ing o 12 o a ian samples and 4 endome ial samples we e included in he mRNA p ofiling analysis. Addi ionally, all 16 sam- ples we e included in he miRNA p ofiling analysis, excep o M2-EA (Table S1). A e he sample collec ion, he issue samples o domes ic sheep we e quickly s o ed in RNAla e eagen (Ambion/Qiagen, Valencia, CA, USA) pe he manu- ac u e ’s ins uc ions and anspo ed o he labo a o y. Fo he Eu opean mouflons, he collec ed samples we e s o ed in RNAla e wi hin 1 h o dea h due o he o es a eas. The issue samples we e hen ans e ed o be s o ed a −80°C in labo a o y un il ex ac ion. We used an RNeasy plus mini ki (Qiagen, Valencia, CA, USA) o ex ac he mRNA and miRNA om he issues acco ding o he manu ac u e ’s p o ocol. The RNA concen a ion and RNA in eg i y numbe we e measu ed using a Bioanalyze 2100 (Agilen Technologies, Waldb onn, Ge many). Lib a ies o mRNAs and miRNAs we e p epa ed using he Illumina’s T uSeq lib a y p epa a ion ki s and sequenced using Illumina Hiseq2000 a he Finnish Mic oa ay and Sequencing Cen e , Tu ku, Finland. The mRNA lib a ies we e sequenced using 100 base-pai (bp) pai ed-end sequencing chemis y, whe eas he miRNAs we e sequenced using he single-end 50 bp app oach. Quali y con ol and mapping Fas QC .0.11.4 (h p://www.bioin o ma ics.bab aham.ac. uk/p ojec s/ as qc/) was used o examine he quali y o he aw RNA-Seq (bo h he mRNA-Seq and miRNA-Seq) eads, and he adap e sequences we e emo ed using he p og am Cu adap .1.9.1 (Ma in 2011). The ob ained immed mRNA-Seq eads we e hen mapped agains he sheep e e ence genome (Oa .4.0) suppo ed by he gene anno a ion file o Ensembl elease 83 using he p og am TopHa .2.0.8b (T apnell e al. 2012). Fo he miRNA-Seq da a, he p og am PRINSEQ-LITE .0.20.4 (Schmiede and Edwa ds 2011) was used o u he miRNA ead p e- p ocessing (including fil e ing o eads con aining ambig- uous eads, da a e o ma ing, and adap e imming). Clean eads be ween 18 and 26 nucleo ides (n ) in leng h we e used o he subsequen analysis. Then, he Bow ie-build ool was used o build an index o he genome alignmen using he sheep genome assembly Oa .3.1 ( p:// p. ensembl.o g/pub/ elease-83/ as a/o is_a ies/dna/), and he clean eads (18–26 n ) we e mapped o he index using he p og am Bow ie .2.2.1 (Langmead e al. 2009). The bioin o ma ics pipeline used o analyze he mRNA-Seq and miRNA-Seq da a is summa ized in Fig. S1. De no o ansc ip ome assembly and anno a ion o Eu opean mouflon We used he T ini y .2.1.1 package (G abhe e al. 2011) o cons uc a de no o ansc ip ome assembly o Eu opean mouflon using he o a ian (M1-OA, M1-OB, M2-OA, M2- OB, M3-OA, and M3-OB) and endome ial (M2-EA and M2-EB) mRNA sequences. The quali y o he de no o assembly was assessed by examining he p opo ion o eads mapped o he assembly. We used he RSEM (RNA- Seq by expec a ion-maximiza ion) ool (Li and Dewey 2011) o es ima e he numbe o RNA-Seq agmen s ha map o each con ig in he indi idual M2 (M2-OA, M2-OB, M2-EA, and M2-EB), and de ec ed he di e en ially exp essed genes be ween he wo ypes o issues, i.e., o a y and endome ium, using he EdgeR package (Robinson e al. 2010). Unigenes ha had |log2( old change)| ≥2 and a alse disco e y a e ≤0.01 we e conside ed significan ly di e en ially exp essed. The di e en ially exp essed uni- genes we e mapped o he NCBI n collec ion da abase using he Blas N (Al schul e al. 1997) ool in he Blas .2.4.0+package (h p://blas .ncbi.nlm.nih.go ) wi h a ma ched s anda d “e- alue =1e-10 and ma ching a io = 100%”. The p og am PANTHER .10 (Mi e al. 2016) was used o pe o m he Gene On ology (GO) and Kyo o Encyclopedia o Genes and Genomes (KEGG) analyses, and he h eshold was se o a P- alue < 0.05 a e he Bon e oni co ec ion. The bioin o ma ics analysis wo k- flow is summa ized in Fig. S2. mRNA di e en ial exp ession and unc ional analysis A e he mapping, he p og ams Cu flinks, Cu me ge, and Cu di in he Cu flinks so wa e package .2.1.1 (T apnell e al. 2012) we e used o assemble he mRNA clean eads Compa a i e mRNA and miRNA exp ession in Eu opean mouflon (O is musimon) and sheep (O is a ies). . . and calcula e he gene exp ession le els based on he agmen s pe kilobase o exon pe million agmen s mapped (FPKM) alues. Fi s , we pe o med he mRNA ansc ip assembly and c ea ed a “me ged ansc ip s.g ” o all o a ian o endome ial samples om one species using he p og ams Cu flinks and Cu me ge. Then, he p og am Cu di was used o calcula e he le el o gene exp ession and pe o m di e en ial exp ession analyses be ween he wo species (i.e., o a y/endome ium: Eu - opean mouflon s. Finnsheep) o wo ypes o issues (i.e., Eu opean mouflon: o a y s. endome ium). The di e en- ially exp essed genes (DEGs) we e defined as genes ha showed |log2( old change)| ≥2 and Padj (i.e. adjus ed P- alue) ≤0.01 in he compa isons. Fu he mo e, we pe - o med GO and KEGG pa hway en ichmen analyses o unc ional genes associa ed wi h he di e en ially exp essed mRNAs using he Da abase o Anno a ion, Visualiza ion, and In eg a ed Disco e y (DAVID) ools .6.8 (Huang e al. 2009). The en iched GO e ms we e conside ed significan a an adjus ed P- alue < 0.05. A wo d cloud was c ea ed o he GO e ms and pa hways using an online Wo dcloud gene a o (h ps://www.jasonda ies.com/wo dcloud/). miRNA quan ifica ion and di e en ial exp ession analysis The miRDeep2 .0.0.7 so wa e package (F iedlände e al. 2012) was used o p edic he pu a i e p ecu so miRNA sequences, and p ecu so miRNAs wi h 10 o ewe eads and miRDeep2 sco es < 5 we e emo ed. Then, he fil e ed miRNAs we e aligned o he miRNA da abase miRBase .21 (h p://www.mi base.o g/) o iden i y ma u e miRNAs ha con ain sheep miRNAs (i.e., miRNAs deposi ed in he miRBase o O. a ies), conse ed miRNAs (i.e., miRNAs ma ching hose in o he mammals in miRBase wi h ewe han ou misma ches), and no el miRNAs ( emaining miRNAs). The di e en ially exp essed o a ian and endome ial miRNAs (including he sheep and conse ed miRNAs) be ween Eu opean mouflon and Finnsheep and be ween he wo ypes o issues (o a y s. endome ium) in Eu opean mouflon we e iden ified using he DESeq .1.22.0 p og am (Ande s and Hube 2010). By applying he h eshold alues o Padj (i.e., adjus ed P- alue) ≤0.05 and |log2( old change)| ≥1, he emaining miRNAs conside ed o be up egula ed o down egula ed in one o he wo species o ypes o issues we e used in he ollowing mRNA–miRNA in eg a ed analyses. miRNA–mRNA egula o y ne wo k analysis The p og am Ta ge Scan .7.0 (Aga wal e al. 2015) was used o p edic he binding o he di e en ially exp essed miRNAs o hei pu a i e a ge s, i.e., he 5′seed egions (2–8 n ) o he miRNAs complemen a y o he conse ed 8- and 7-me si es in he 3′-un ansla ed egion o he mRNA (Lewis e al. 2005). miRNAs ypically unc ion in pos - ansc ip ional egula ion by supp essing o silencing spe- cific a ge genes. Thus, we sea ched o in e se co ela ions be ween he exp ession le els o he miRNAs and hei p edic ed a ge genes, i.e., be ween he up egula ed miR- NAs and down egula ed mRNAs, and ice e sa. To na ow he lis o po en ial miRNA– a ge in e ac ions, he pu a i e a ge genes we e o e lapped wi h he DEGs in he mRNA p ofiling da a. Then, a ge genes wi h a con ex ++ sco e pe cen ile lowe han 50 we e emo ed. Fu he mo e, he up egula ed and down egula ed miRNAs and hei a ge genes we e included in he ne wo k isualiza ion analysis using he p og am Cy oscape .3.4.0 (Shannon e al. 2003). Only he op 10 down egula ed miRNAs we e conside ed co e miRNAs and used o build he main egula o y ne - wo k. Fu he mo e, GO and pa hway en ichmen analyses o he miRNA a ge genes we e conduc ed using DAVID ools .6.8 (Huang e al. 2009). GO e ms and pa hways wi h an adjus ed P- alue < 0.05 we e conside ed sig- nifican ly en iched. A flowcha o he mRNA–miRNA in eg a ed analysis is p o ided in Fig. S1. To confi m he c edibili y o he in e ed miRNA and a ge gene pai s, we sea ched he miRTa Base (h p://mi a base.mbc.nc u.edu. w/; Chou e al. 2016) da abase and p e ious s udies o di e en ially exp essed miRNAs, p e- dic ed a ge genes, and miRNA–gene pai s consis en wi h hose iden ified in his s udy. Resul s and Discussion mRNA and miRNA p ofiling In o al, we gene a ed app oxima ely 127.1 Gb mRNA and 21.3 Gb miRNA aw da a. A e fil e ing, we ob ained a la ge numbe o high-quali y eads: 1083.6 million mRNA and 56.8 million miRNA eads in Eu opean mouflon, and 934.6 million mRNA and 127.8 million miRNA eads in Finnsheep. Mos high-quali y mRNA and miRNAs ha e been mapped o he O.a ies e e ence genome Oa _ .4.0 (Table S2). Rega ding he mRNA exp ession, mo e han 89% o he genes we e exp essed a <50 FPKM and only <1.1% o he genes we e exp essed a >500 FPKM in all ou issue-species combina ions (Table S3). Al oge he , 23 557 common genes we e exp essed in he o a ies and 19 208 common genes we e exp essed in he endome ia (Table S3). In he de no o assembly o he Eu opean mouflon an- sc ip ome, 1083.6 million clean mRNA eads we e assem- bled in o 844 129 unigenes, wi h an N50 leng h o 700 bp Ji Yang e al. and an a e age leng h o 562 bp (Table 1). The obse ed leng h o N50 is compa able o ha epo ed in domes ic sheep (e.g., 508–1482 bp) (Yue e al. 2015; Zhang e al. 2015), hus highligh ing he high quali y o he assembled Eu opean mouflon ansc ip ome. Based on he an- sc ip ome, we iden ified 5739, 5831, 5897, and 5981 anno a ed DEGs in he ou compa isons be ween he o a y and endome ium issues (M2-OA s. M2-EA, M2-OA s. M2-EB, M2-OB s. M2-EA, and M2-OB s. M2-EB) in Eu opean mouflon, espec i ely. Howe e , when he ana- lyses we e pe o med using he sheep e e ence genome (Oa .4.0), we ound a ewe DEGs (e.g., 1925 DEGs o M2-O-B s. M2-E-B). Despi e he possible implica ion o using he O. a ies genome as a e e ence, he obse a ion sugges s ha he de no o app oach is accu a e and ha he quali y o he de no o assembly is good. We iden ified 373 ma u e miRNAs in he o a ies and 192 ma u e miRNAs in he endome ia, including 89 and 59 sheep miRNAs, 174 and 104 conse ed miRNAs, and 110 and 29 no el miRNAs in he wo issues, espec i ely (Table S4). High le els o exp ession we e obse ed in he sheep miRNAs, ollowed by conse ed miRNAs and no el miRNAs (Table S4). O he 122 no el miRNAs (Table S5), 20 miRNAs we e also iden ified as no el miRNAs and alida ed in ou p e ious s udy (Hu e al. 2016). The ch omosomal dis ibu ion o he sheep and conse ed miRNAs indica ed he p esence o la ge numbe s o miR- NAs on ch omosomes 18 and X, whe eas none o he miRNAs we e exp essed on ch omosome 8 based on he o a ian da a and ch omosomes 8, 17, and 25 based on he endome ium da a (Table S6). mRNAs and miRNAs showing he g ea es exp ession O he genes exhibi ing he g ea es exp ession (FPKM > 3000) in o a y, 7 genes in he Finnsheep and 15 genes in he Eu opean mouflon encode ibosomal p o ein (Fig. 1a and Table S7). In he endome ium, only 1 op gene iden- ified in he Finnsheep and 22 op genes iden ified in he Eu opean mouflon encode ibosomal p o ein in ol ed in he ibosome pa hway (Fig. 1b and Table S7). The ibosome pa hway has been epo ed o be among he mos significan pa hways in ol ed in ea ly emb yogenesis (Be es e al. 2017). The obse ed la ge numbe o ibosome amily genes in bo h issues o Eu opean mouflon could be one explana ion o he low ep oduc i e success o he wild sheep (i.e., seasonal es us, 1.5 yea s o age o sexual ma u a ion and one lamb pe li e ; Cugnasse e al. 1985; Ga el e al. 2005) compa ed o he p olificacy o domes ic sheep, including Finnsheep (i.e., ou -o -season b eeding, 6 mon hs o age o sexual ma u a ion and h ee lambs pe li e on a e age; Maijala 1996; Li e al. 2011; Mullen and Han ahan 2014). Mo eo e , he obse ed smalle numbe o ibosome amily genes exp essed a a highe le el in Finnsheep han in Eu opean mouflon is likely due o he e olu ion o ep oduc i e ai s be ween he wo species du ing he domes ica ion and long- e m a ificial selec ion p ocesses. Addi ionally, his obse a ion could be ela ed o he en i onmen o he mouflon being ha she han he domes ic en i onmen o he Finnsheep; in such a ha sh en i onmen , i may be beneficial o ha e ewe o sp ing. O he 10 mos highly exp essed o a ian miRNAs, 7 (miR-10b,miR-143,miR-26a,miR-27b,le -7 ,miR-191, and miR-22-3p) miRNAs a e common be ween Eu opean mouflon and Finnsheep (Table S8). O he endome ial miRNAs, h ee (miR-10b,miR-100, and miR-34c-5p) miR- NAs a e highly exp essed in bo h species (Table S8). We also obse ed issue- and species-specific miRNAs (Table S9), which could possibly be associa ed wi h pa icula unc ions. Fo ins ance, miR-181a, which is specifically and abundan ly exp essed in he Finnsheep o a y, was demon- s a ed o be associa ed wi h p egnancy complica ions and was mo e highly exp essed in no mal e m han in p e- e m deli e y (Mayo -Lynn e al. 2011). Addi ionally, miR-186, which is specifically and abundan ly exp essed in he Eu - opean mouflon endome ium, was obse ed o be mo e highly exp essed in endome ial cance (Zhou e al. 2008). Di e en ially exp essed mRNAs and miRNAs In o al, 192 and 1065 DEGs esponsible o o a ian and endome ial mRNAs we e iden ified be ween he wo spe- cies, espec i ely. O hese DEGs, 149 o a ian and 492 endome ial DEGs we e up egula ed in Eu opean mouflon, Table 1 Summa y o cons i uen da a o T ini y-assembled Eu opean mouflon ansc ip ome Desc ip ion Numbe Be o e imming Raw eads 1 083 877 020 A e imming Clean eads 1 083 600 777 A e age leng h o clean eads (bp) 101 N pe cen age 0 GC con en (%) 48.92 A e assembly To al T ini y “genes”844 129 To al T ini y ansc ip s 1 102 841 S a is ics based on all ansc ip con igs Mean leng h o con igs (bp) 947.39 N50 (bp) o con igs 2169 S a is ics based on unigenes Mean leng h o unigenes (bp) 562.14 N50 (bp) o unigenes 700 Compa a i e mRNA and miRNA exp ession in Eu opean mouflon (O is musimon) and sheep (O is a ies). . . while 43 o a ian and 573 endome ial DEGs we e up e- gula ed in Finnsheep (Tables S10 and S11). Addi ionally, we iden ified 53 and 5 di e en ially exp essed o a ian and endome ial miRNAs be ween he wo species, espec i ely (Table 2). Based on he di e en ially exp essed mRNAs and miRNAs in he issues, hea maps we e cons uc ed (Fig. 2) and showed dis inc exp ession pa e ns be ween Eu opean mouflon and Finnsheep, indica ing ha he e o- lu ion o he gene- and pos - ansc ip ional egula ion could ha e con ibu ed o he imp o ed ep oduc i e pe o mance du ing he domes ica ion p ocess. O he up egula ed genes in he o a y and endome ium issues om he Finnsheep, 4 genes (CYP19,PTGS2,FST, and PAPSS2) a e associa ed wi h ea ly ma u i y, 14 genes (EREG,INHBA,SERPINA5,SPP1,LIPG,PTGS2,AREG, TDGF1,FST,C4BPA,SERPINE2,TRPV3,LDLR, and TINAGL1) a e associa ed wi h ecundi y, and wo genes (CYP19 and SERPINE2) a e associa ed wi h non-seasonal b eeding. In addi ion, we compa ed he 192 o a ian DEGs wi h he QTLs o he ep oduc i e ai s in he Sheep QTL da abase (h ps://www.animalgenome.o g/cgi-bin/QTLdb/ OA/), and ound 8 genes (NME9,LRRIQ4,SPATC1L, PRDM15,PEX5L,COL6A6,RARRES1, and SPOCK1; Table S12) o e lapped wi h p e iously iden ified QTLs associa ed wi h ai s such as ep oduc i e seasonali y and o al lambs bo n. Simila ly, we ound 73 endome ial DEGs ha o e lapped wi h he epo ed QTLs o he ep oduc i e seasonali y, and 1 gene, i.e., RASL11B, was o e lapped wi h he QTL ch omosome 6 as ollows: 68295117–68445447 o he o al lambs bo n ai (Table S13). GO and KEGG en ichmen s o di e en ially exp essed mRNAs In he unc ional en ichmen analysis o he 192 o a ian DEGs be ween Eu opean mouflon and Finnsheep, he op en iched GO ca ego ies we e mainly ele an o ep oduc- ion, including e iliza ion, ep oduc i e de elopmen al p ocess, ep oduc i e p ocess, oogenesis, and o a ian ol- licle de elopmen (Fig. S3a and Table S14). Fe iliza ion e e s o he usion o game es ha ini ia es he de elopmen o a new indi idual o ganism (Me iam-Webs e 2017). Du ing he p ocess, he ZP p o eins (ZP1 o ZP4) in he zona pellucida unc ion o su ound he egg, bind he spe m and media e species-selec i e spe m–oocy e in e ac ions (Suzuki e al. 2015; A ella e al. 2016). ZP2 pep ide beads ha e been shown o ma kedly inhibi he e iliza ion o o ula ed eggs, p e en spe m binding in mice and humans (Gup a e al. 2012; A ella e al. 2016) and induce in e ili y in emale mice (A ella e al. 2016). ASTL, which is in ol ed in he nega i e egula ion o e iliza ion by he ZP2 and ZP4 genes, plays a defini i e ole in ensu ing monospe mic e iliza ion by encoding o as acin, which is a co ical g anule p o ease (Bu ka e al. 2012). Thus, he significan ly up egula ed exp ession o he ASTL,ZP2,ZP3, and ZP4 genes in he Eu opean mouflon o a y (Table S10) may indica e ha hese genes play a ole in he low e ili y o wild sheep compa ed o he high p olificacy o domes ic sheep (Fig. 3). The GO e ms o he ep oduc i e de elopmen p ocess and o a ian ollicle de elopmen included he genes EREG, INHBA,SERPINA5, and SPP1, which a e significan ly up egula ed in Finnsheep (Tables S10 and S14). EREG is a membe o he epide mal g ow h ac o amily, which plays an impo an ole in media ing he lu einizing ho mone- induced o ula ion p ocess in p eo ula o y ollicles (Ams e dam 2010; Kim e al. 2011). INHBA plays an essen ial ole in ea ly oocy e and ollicula de elopmen (Thomas e al. 2003) and has a significan e ec on he li e size in sheep (Leyhe e al. 1994) and goa (Zi e al. 2012). SERPINA5, which encodes a se ine p o ease inhibi o , may be in ol ed in ollicula g ow h and a esia in he bo ine Fig. 1 Numbe o genes exhibi ing he g ea es exp ession (FPKM > 3000) in he Finnsheep and Eu opean mouflon o a ies (a) and endome ia (b), and hei o e lap wi h genes in he ibosome pa hway Ji Yang e al. o a ian ollicle (Hayashi e al. 2011). SPP1 is loca ed on sheep ch omosome 6 close o he high p olificacy gene FecB in sheep (Mon gome y e al. 1993; Mon gome y e al. 1994). Addi ionally, he genes AMH,TDRD5,SOHLH1, and NOBOX a e up egula ed in Eu opean mouflon (Tables S10 and S14). AMH is a cen al influencing ac o in gonad de elopmen , and could inhibi ge m cell p oli e a ion and di e en ia ion in bo h sexes (Du linge e al. 2002). Recen clinical da a ha e indica ed ha se um AMH le el is a use ul indica o o o a ian aging and p ema u e o a ian ailu e in women (de Ve e al. 2002; Alipou e al. 2015). The TDRD5 p o ein is an e olu iona ily conse ed p o ein in ol ed in ge m cell de elopmen . Injec ion o TDRD5- deficien ound spe ma ids in o oocy es esul in e ile o sp ing (Yabu a e al. 2011). SOHLH1 is a c i ical eg- ula o o oogenesis and can down egula e NOBOX, which unc ions o dis up olliculogenesis (Pangas e al. 2006). The e o e, ou esul s sugges ha he up egula ed Table 2 Lis o significan ly di e en ially exp essed miRNAs be ween Finnsheep and Eu opean mouflon issues (o a y and endome ium) miRNA log2( old change)aPadjbmiRNA log2( old change)aPadjb O a y (Eu opean mouflon up o Finnsheep down) O a y (Eu opean mouflon down o 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-30 1.65 5.96E-05 miR-10b 1.81 8.11E-03 le -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 le -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 Endome ium (Eu opean mouflon up o Finnsheep down) Endome ium (Eu opean mouflon down o 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 h eshold o old change is se o log2( old change) > 1 o <−1 bThe s a is ical significance is e alua ed as Padj (i.e., adjus ed P- alue) < 0.05 Compa a i e mRNA and miRNA exp ession in Eu opean mouflon (O is musimon) and sheep (O is a ies). . . exp ession o hese genes (i.e., EREG,INHBA,SERPINA5, and SPP1 in Finnsheep, and AMH,SOHLH1,TDRD5, and NOBOX in Eu opean mouflon) is mos likely associa ed wi h he hype -p olificacy o domes ic sheep (e.g., Finn- sheep) and he single lambing o wild sheep (e.g., Eu opean mouflon; Fig. 3). We also obse ed o he op GO ca ego ies Fig. 2 Hea maps o di e en ially exp essed mRNAs and miRNAs be ween Finnsheep and Eu opean mouflon. aHea map showing he exp ession in ensi y o 1065 di e en ially exp essed genes (DEGs) based on a compa ison o endome ia om he wo species. bHea map showing he exp ession in ensi y o 192 DEGs based on a compa ison o o a ies om he wo species. c Hea map showing he exp ession in ensi y o 5 di e en ially exp essed miRNAs based on a compa ison o endome ia om he wo species. dHea map showing he exp ession in ensi y o 53 di e en ially exp essed miRNAs based on a compa ison o o a ies om he wo species. Hea maps we e gene a ed using he R hea map package .1.0.8 (h p://c an. -p ojec .o g/web/pa ckages/phea map/), which is based on he k-means clus e ing algo i hm and he Boolean alues o de e mine i he ows (i.e., mRNAs o miRNAs) should be clus e ed. The impo an genes associa ed wi h he ep oduc i e ai s such as ea ly ma u e, ecundi y, non- seasonal b eeding, and o al lambs bo n a e labeled in he figu e Ji Yang e al. ela ed o immune p ocesses such as he humo al immune esponse and humo al immune esponse media ed by ci - cula ing immunoglobulin (Table S14). The KEGG analysis o he 192 o a ian DEGs iden ified h ee significan (adjus ed P- alue < 0.05) pa hways, i.e., o a ian s e - oidogenesis, complemen and coagula ion cascades, and MAPK signaling pa hways (Table S14). The MAPK (i.e., p38αMAPK) pa hway has been shown o play an essen ial ole in emale ep oduc ion, pa icula ly in ollicula de elopmen (Hu e al. 2017) (Fig. 3). In he unc ional en ichmen o he endome ial DEGs be ween he wild and domes ic sheep, he GO e ms o he 492 up egula ed genes (RARRES2,LTBP1, and TGFB3)in Eu opean mouflon we e mos significan ly en iched in he de elopmen al p ocess, such as ana omical s uc u e de elopmen (Fig. S3b and Table 3), whe eas he op en iched GO e ms o he 573 up egula ed genes (PTGS2, LDLR, and SREBF1) in Finnsheep included biosyn he ic and me abolic p ocesses, such as lipid me abolic p ocess (Fig. S3b and Table 3). An ea ly s udy in es iga ing he bo ine endome ial ansc ip ome indica ed an ea lie shi om p oli e a ion o me abolism a ea ly dies us in he cow (Mesqui a e al. 2015). Cell p oli e a ion was en iched in he small ollicle endome ium, whe eas biosyn he ic and me abolic p ocesses we e en iched in he la ge ollicle endome ium, which p esen ed highe p oli e a i e ac i i y in he luminal epi helium, glandula epi helium, and s oma han in he small ollicle endome ium (Mesqui a e al. 2015). Addi ionally, he la ge ollicle endome ium is an ac i e endome ial pheno ype ha is associa ed wi h he op imal u e ine en i onmen (Mille and Moo e 1976; Mesqui a e al. 2015). Thus, Finnsheep could ha e a la ge ollicle endome ium and be e ecep i i y o he u e us han Eu opean mouflon (Fig. 3). In he KEGG analysis o he 1065 endome ial DEGs, he op en iched pa hways we e en iched in he ex acellula ma ix (ECM)– ecep o in e ac ion, s e oid biosyn hesis, and biosyn he ic and me abolic p ocesses (Table S15). The ECM, which is a complex ma ix o biological mac omolecules, can ans e signals o cells ia su ace ecep o s du ing cell adhesion and o he a ious signaling ansduc ion pa hways. Addi ionally, he ECM can send signals o he cy oplasm and nucleus o influence gene exp ession o cellula ac i i ies (Alldinge e al. 2006). Thus, adhesion and connec ion mechanisms could play impo an oles in he endome ium a he molecula le el, and di e en cell adhesion may occu be ween wild and domes ic sheep. In addi ion, he di e en ially exp essed mRNAs be ween he o a ian and endome ial issues om he Eu opean mouflon and hei unc ional analysis based on he e e ence genome and de no o app oaches a e p esen ed in he Supplemen a y Resul s and Tables S16–S18. Pos - ansc ip ional egula o y ne wo k o miRNAs and a ge genes We ob ained 100, 64, 17, and 51 a ge genes o he up- and down egula ed miRNAs in pa icula issues (Table 4). In he main down egula ed ne wo k o he Finnsheep o a y, which comp ised 10 down egula ed miRNAs and 57 up e- gula ed a ge genes (Fig. 4a), se e al impo an miRNAs Fig. 3 Func ional oles o he main mRNAs and miRNAs in he endome ium and o a y o O is species. O a ian and endome ial ne wo k o di e en ially exp essed mRNAs and op 10 Finnsheep down egula ed (i.e., Eu opean mouflon up egula ed) miRNAs and a ge genes a e shown, along wi h hei unc ional oles associa ed wi h ep oduc ion (bold ype). Regula ion and al e a ion o hese unc ional oles a e esponsible o he e olu ion o high p olificacy ai s. Red and blue a ows indica e down egula ion and up egula ion o he ep oduc i e o gans, espec i ely Compa a i e mRNA and miRNA exp ession in Eu opean mouflon (O is musimon) and sheep (O is a ies). . .