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). . .