A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 1 Transcriptome profiling of Finnsheep ovaries during out-of-season breeding period Kisun Pokharel1,2, Jaana Peippo1, Göran Andersson3, Meng-Hua Li4 and Juha Kantanen1,2 1Green Technology, Natural Resources Institute Finland (Luke), Myllytie 1, FI-31600 Jokioinen, Finland 2Department of Biology, University of Eastern Finland, Yliopistonranta 1 E, FI-70211 Kuopio, Finland 3Department of Animal Breeding and Genetics, Swedish University of Agricultural Sciences, Gerda Ulls Väg 26, 75007 Uppsala, Sweden 4CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences (CAS), Beichen West Road No. 1–5, Chaoyang District, 100101 Beijing, China email: menghua.li@ioz.ac.cn Finnsheep is one of the most prolific sheep breeds in the world. We sequenced RNA-Seq libraries from the ova- ries of Finnsheep ewes collected during out of season breeding period at about 30X sequence coverage. A total of 86966348 and 105587994 reads from two samples were mapped against latest available ovine reference genome (Oarv3.1). The transcriptome assembly revealed 14870 known ovine genes, including the 15 candidate genes for fertility and out-of-season breeding. In this study we successfully used our bioinformatics pipeline to assemble the first ovarian transcriptome of Finnsheep. Key words: bioinformatics pipeline, mRNA, prolific sheep breed Introduction The native Finnsheep breed demonstrates several economically important and physiologically interesting traits such as large litter size (i.e. multiple births) and out-of-season breeding and, thus, has been exported to more than 40 countries (Maijala 1988). The large litter size in Finnsheep (on average 2.4–2.9 lambs) is most likely as- sociated with increased ovulation rate during the oestrus cycle and successful subsequent embryonic and foetal development (Maijala and Österberg 1977, Fahmy and Dufour 1988, Li et al. 2011, Vinet et al. 2012). Mutations in functional genes such as bone morphogenetic protein 15 (e.g. BMP15 on OARX), bone morphogenetic protein receptor, type 1B (BMPR1B on OAR6) and growth differentiation factor 9 (GDF9 on OAR5) have been shown to account for the increase in ovulation rate and litter size in sheep (McNatty et al. 2005, Vinet et al. 2012, Mullen et al. 2013, Våge et al. 2013). The majority of Finnsheep ewes also have an exceptional ability to breed out-of-season, i.e. they can conceive all the year around instead of only once in autumn, typical for most breeds, allowing year-round breeding (Sormunen- Cristian and Suvela 1999). Functional genes accounting for the off-season oestrus in Finnsheep have not yet been identified although a few candidate genes for the trait have been found in some other prolific sheep breeds. For example, Chen et al. (2012) identified a set of candidate genes for out-of-season breeding in the Hu sheep (an in- digenous Chinese breed), such as neurotrophic tyrosine receptor type 2 (NTRK2), nidogen 1 (NID1) and serine pro- tease inhibitor clade E member 2 (SERPINE2). Moreover, Notter et al. (2003) and Mateescu and Thonney (2010) found another candidate gene, melatonin 1A receptor (MTNR1A), for out-of-season breeding in the Dorset sheep. Earlier investigations showed that functional genomics elements underlying complex traits can be efficiently ex- plored by whole-transcriptome analysis, using next-generation sequencing technologies (Twine et al. 2011, Huss- eneder et al. 2012, Ramayo-Caldas et al. 2012, Xu et al. 2013). We have commenced a sheep fertility study which is based on whole-genome transcriptome and single nucleotide polymorphism (SNP) marker analyses. In the pre- sent study, we present our first results on transcriptome profiling of Finnsheep ovaries and test the applicability of experimental design and methods prior to analyses of a larger, main data set. The ovaries were collected dur- ing the out-of-season breeding period (June) from two Finnsheep ewes. To the best of our knowledge, this is the first transcriptome profiling study on the exceptionally prolific Finnsheep breed, which represents a valuable ani- mal genetic resource for the global sheep industry. Manuscript received August 2014 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 2 Materials and methods Sample collection Animal handling procedures and sample collection were conducted following both the Finnish laws and EU di- rectives regarding animal experimentation and occupational health. Ovarian mRNA expression profiles were in- vestigated in two Finnsheep ewes. The two ewes were born in January 2011 on the same farm located in west- ern Finland and were closely related (coefficient of relation between the ewes was 0.375 as estimated from the pedigree records). One ovary from each of the ewes was collected after slaughter during June 2012. The samples were decoded and hereafter are referred as S48 and S50. RNA isolation and sequencing Ovaries were dissected into 4 sections and stored in RNAlater reagent (Qiagen, Valencia, CA, USA) at –20 °C until use. Total RNA was extracted using the RNeasy plus mini kit (Qiagen, Valencia, CA, USA). RNA concentration and RIN (RNA Integrity Number) values were checked using a Bioanalyzer 2100 (Agilent Technologies). cDNA libraries were prepared and sequenced using Illumina HiSeq 2000 (Illumina, San Diego, CA) at the Institute of Molecular Medicine Finland as detailed in Koskela et al. (2012), generating 100-bp paired-end reads. The RNA-Seq sequenc- ing data have been deposited in the Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra/) with the accession number SRX791949. Transcriptome assembly and annotation The used bioinformatics pipeline is illustrated in Figure 1. The quality of reads was assessed using FastQC version 0.10.1 (Simon Andrews). The sheep reference genome (Oarv 3.1) and reference genes were downloaded from En- sembl database (Ensembl Release 74). We mapped the obtained reads to the ovine reference genome and gene sequences using Tophat v2.0.9 (Kim et al. 2013). Cufflinks v2.1.1 (Trapnell et al. 2010) was used to assemble tran- scripts from the RNA-Seq alignments. In addition, two different assemblies for each sample were compared and combined with reference transcriptome to identify unknown transcripts in the data using the Cuffcompare pro- gram of Cufflinks. All reads that exclusively mapped to ovine genes were quantified using HTSeq-count software package (Anders et al. 2015). Gene expression levels of read counts obtained from HTSeq-count were measured in the R statistical environment using the DESeq Bioconductor package (Anders and Huber 2010). Fig. 1. RNA-Seq pipeline. Flowchart illustrates the bioinformatics pipeline used to analyse the RNA-Seq data. Arrows indicate data flow; analysis steps are marked by hexagons; rectangles indicate output data and rounded rectangles indicate tools and programs. A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 3 We used Blast2GO (Conesa et al. 2005) for retrieving Kyoto Encyclopedia of Genes and Genomes (KEGG) path- ways and analyzing annotation data. Ensembl BioMart (Kinsella et al. 2011) was used to retrieve Gene Ontology (GO) terms associated with all expressed genes and updated annotation for Ensembl IDs with missing informa- tion. We annotated unknown transcripts, i.e. all transcripts that did not match to sheep genes, by implementing BLAST search and InterPro query in Blast2GO. Genes associated with fertility and out-of-season breeding were analyzed in detail. Both the nucleotide and amino acid variants in three major candidate genes (GDF9, BMP15 and BMPR1B) were identified and assessed visually using the program Integrative Genomics Viewer (IGV) version 2.3.26 (Thorvaldsdóttir et al. 2013). Results Transcriptome sequencing and assembly A total of 192554342 RNA-Seq reads representing the two studied samples, S48 (86966348 reads) and S50 (105587994 reads), were generated by Illumina sequencing. The overall quality of the reads was very good with an average quality score across all bases greater than 28. Thus, we did not perform any filtration and trimming of the data. GC contents for S48 and S50 were 54% and 55%, respectively. High quality reads were aligned to the ovine genome using Tophat v2.0.9. Without any filtering for base quality scores, 68.6% and 68.2% of the RNA-seq reads from the S50 and S48, respectively, were mapped to the ovine reference genome. Gene expression profiling The recent ovine reference genome (Oarv3.1, December 2013) of Ensembl was used for annotating and counting reads that mapped to genomic features. Altogether 14875 genes were expressed in both samples with the base mean expression values of >10. The ten most highly expressed genes (COL1A1, COL1A2, COX1, ENSOARG00000000004, ENSOARG00000006149, COL3A1, EEF2, ACTB, INHA and SERPINE2) contributed to 9.1% (2272788) of the reads in the S48 and 11.5% (3691546) of the reads in the S50, respectively. INHA, a candidate gene for fertility and SERPINE2, a proposed candidate gene for out-of-season breeding, were the two exceptional highly expressed genes in S50. Out of 14875 Ensembl genes expressed in this study 2261 did not have genes names. We converted IDs with miss- ing gene names into RefSeq protein ID using BioMart and extracted the corresponding gene names from the Ref- Seq database. We were able to find gene names for 1465 Ensembl IDs including 562 genes with uncertain func- tions (those starting with LOC and 202 duplicated IDs). Gene annotation Genes expressed in the data were associated with 80261 GO terms which can be classified into three categories biological processes, cellular component and molecular functions (Fig. 2). 52 % 25 % 23 % Biological processes Cellular components Molecular function Fig. 2. Gene Ontology categories. The figure represents the pie chart distribution of three GO categories, namely, biological processes, cellular components and molecular function. All 14875 genes expressed in the data were used for the ontology analysis. A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 4 The majority of these terms (41663) were associated with biological processes whereas 19895 terms were related to cellular components and 18703 to molecular function. We further expanded three main categories into lower level classification to understand what kind of processes and functions were linked with the genes expressed in the data. From the list of top ten terms in each category, it was found that protein binding is the most common ontology term which is represented by 5572 genes (Table 1). The analysis revealed that expressed genes were mainly involved in differentially regulating transcription processes followed by developmental processes such as cell proliferation and embryonic development. Similarly, nucleus is the most frequent cellular component term annotated to expressed genes. Table 1. Gene Ontology analysis. The table represents top ten GO terms for biological processes, cellular component and molecular function. No GO Terms GO ID No. of Genes Biological Process 1 positive regulation of transcription from RNA polymerase II promoter GO:0045944 485 2 negative regulation of transcription from RNA polymerase II promoter GO:0000122 357 3 positive regulation of transcription, DNA-dependent GO:0045893 349 4 negative regulation of transcription, DNA-dependent GO:0045892 297 5 protein phosphorylation GO:006468 240 6 regulation of transcription, DNA-dependent GO:0006355 222 7 negative regulation of cell proliferation GO:0008285 216 8 negative regulation of apoptotic process GO:0043066 206 9 in utero embryonic development GO:0001701 205 10 positive regulation of cell proliferation GO:0008284 200 Cellular Component 1 nucleus GO:0005634 2442 2 cytoplasm GO:0005737 2375 3 mitochondrion GO:0005739 1048 4 plasma membrane GO:0005886 938 5 cytosol GO:0005829 455 6 membrane GO:0016020 437 7 golgi apparatus GO:0005794 394 8 nucleolus GO:0005730 369 9 endoplasmic reticulum GO:0005783 346 10 centrosome GO:0005813 289 Molecular Function 1 protein binding GO:0005515 5572 2 DNA binding GO:0003677 376 3 protein homodimerization activity GO:0042803 351 4 identical protein binding GO:0042802 344 5 sequence-specific DNA binding transcription factor activity GO:0003700 242 6 chromatin binding GO:0003682 221 7 protein kinase binding GO:0019901 191 8 transcription factor binding GO:0008134 163 9 protein heterodimerization activity GO:0046982 162 10 sequence-specific DNA binding GO:0043565 160 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 5 The obtained transcripts and genes were related to 129 KEGG pathways. Table 2 shows the ten most abundant pathways with corresponding number of sequences and enzymes. The results indicated that the purine metabo- lism pathway comprised the highest number of sequences (566 in the S48; 531 in the S50) and enzymes (58 in the S48; 54 in the S50) in both samples. The number of enzymes in most of the pathways was same although the number of sequences was different. Altogether 5606 unknown transcripts were present in the data with 3060 sequences known in other species and the rest 2546 did not match any sequences in the non-redundant (nr) sequence database. Out of 3060 transcripts that shared sequence similarity with genes found in other species, gene ontology terms were available for 374 transcripts. Interestingly, BLAST searching showed that as many as 225 transcripts were part of the Orf virus (ORFV) genome (Fig. 2). In the list of transcripts, the ORFV was the third most common genome after bovine and yak. Pres- ence of the viral genes was verified using BLAST search for each assembled transcriptome (S48 and S50) separately. Fig. 3. Distribution of unknown transcripts that did not align to the ovine genes in Ensembl database. Table 2. KEGG pathway analysis. The table shows top ten KEGG pathways and corresponding number of gene sequences (#Seqs) and enzymes (#Enzs) in the data. Pathways S48 S50 # Seqs # Enzs # Seqs # Enzs Purine metabolism 566 58 531 54 Arginine and proline metabolism 420 31 376 31 beta-Alanine metabolism 385 17 331 17 Pyrimidine metabolism 234 37 227 38 Phosphatidylinositol signalling system 216 22 207 22 Glycerophospholipid metabolism 170 30 172 30 Inositol phosphate metabolism 165 25 152 25 Glycerolipid metabolism 114 16 119 16 Lysine degradation 110 17 104 18 Sphingolipid metabolism 109 21 99 21 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 6 Analysis of candidate genes We further examined the expression level of the 17 putative candidate genes, ten for fertility and seven for out-of-season breeding (Table 3). Genes such as INHA, INHBA and SERPINE2 have differential mean val- ues and seem to be up-regulated in S50. However, we were unable to generate any conclusions regard- ing these genes due to the limited number of samples. B4GALNT2 one of the candidate genes for high litter size, and SERPINE2, another candidate gene for out-of-season breeding were not expressed in the data. Mu- tations in GDF9, B4GALNT2, BMPR1B and BMP15, have been shown to affect litter size and species-specif- ic rates of ovulation (Davis 2005, Fogarty 2009, Demars et al. 2013, Våge et al. 2013). Compared to the ovine reference genome, five SNP mutations were present in GDF9. At the amino acid level, instead of the V371M mutation recently reported by Våge et al. (2013), we found the K241E substitution. Similarly, five SNPs were detected in the BMPR1B gene with three non-synonymous mutations (Q249R, K299E and E461G). Howev- er, the BMP15 gene was found to be monomorphic and no causative mutations were detected in the data. Discussion The ovary is one of the key tissues in a complex regulatory network affecting fertility and breeding in mammalian species. We profiled ovarian transcriptomes using biopsies of two ewes of the native Finnsheep, one of the most fertile sheep breeds in the world. Using the recent Ensembl annotation, we identified a total of 14875 genes ex- pressed in the data. The genes expressed in the data represented the majority of the genes found in ovaries of other mammals such as human and mouse. This result is in good agreement with the recent study by Bonnet et al. (2013), where 15349 genes were identified in the sheep juvenile ovary by applying a laser capture dissection approach. Due to smaller sample size, although we found genes which were differentially expressed between two samples, no further analysis were done. However, these findings will be considered in future experiments. The samples were collected during the out-of-season breeding period. All six candidate genes - NTRK2, NID1, FOXO1, HTRA1, SERPINE2, PPAP2B - which were putatively associated with out-of-season breeding capability (Chen et al. 2012) were expressed in the samples. These putative genes associated with out-of-season breed- ing have not yet been well supported by experimental validations and thus the results need to be supplemented with more samples and experimental studies. In addition to the genes associated with out-of-season breeding, we were interested in the expression patterns of fertility genes and their polymorphisms in the Finnsheep. Nine of the ten putative genes (GDF9, BMP15, BMPR1B, INHA, INHBA, PRLR, FSHR, ESR1, ESR2) associated with fertil- ity were expressed in both ewes, but the recently suggested candidate gene B4GALNT2 was not expressed. In the two samples, the causative mutations at the fertility genes GDF9, BMP15 and BMPR1B affecting ovulation rate Table 3. Expression levels of the candidate genes. The first nine genes are associated with fertility and the rest with out- of-season breeding. Base mean values are the weighted mean values computed by DESeq and p adj is the adjusted p value. Genes Base Mean values (S48) Base Mean values (S50) p adj GDF9 295.272 256.862 0.870 BMP15 25.041 41.213 1.646 BMPR1B 4695.13 6711.946 1.429 INHA 52004.372 544759.278 10.475 INHBA 5912.7 103335.0 17.476 FSHR 1117.4 3537.6 3.166 ESR1 1487.8 1417.5 0.953 ESR2 430.9 588.5 1.365 PRLR 73.035 88.176 1.207 NTRK2 2600.062 1183.672 0.455 NID1 21228.276 15125.120 0.712 FOXO1 7541.432 27351.922 3.627 HTRA1 3034.101 3641.108 1.200 SERPINE2 26099.740 169539.172 6.496 PPAP2B 11195.291 11909.559 1.064 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 7 and litter size (Vinet et al. 2012, Våge et al. 2013) were not found. This observation needs to be verified with a larger dataset, but the results support the hypothesis that the mutations at major genes may not explain the ex- ceptional prolificacy of the Finnsheep. In addition to the candidate genes, the data provide insights into biological pathways in the sheep ovary during the out-of-season breeding period. Follicles, the functional units of the ovary, are made up of an oocyte surrounded by one or more layers of somatic cells, namely theca cells and granulosa cells. Early stage development of follicles is mediated by two members of the transforming growth factor β (TGF-β) family, GDF9 and BMP15 (Dong et al. 1996, Carabatsos et al. 1998, Galloway et al. 2000). Follicle development progresses through successive stages via bi-directional interaction between oocyte-granulosa cells and granulosa cells – theca cells (Eppig 2001), featuring again the members of TGF-β family. Several members of the TGF-β were expressed in the data. These include the prototypic TGF-β subfamily (TGF-β3, TGF-β1, TGF-β1I1, TGF-β1R1, TGF-β2, TGF-βR2, TGF-βI, TGF-βR3), an extensive BMP subfamily (BMP15, BMP1, BMPR1B, BMP6, BMPR2, BMP7, BMP2, BMP3, BMP4), the GDF subfamily (GDF11, GDF6, GDF9, GDF10, GDF3, GDF1), the inhibin subfam- ily (INHA, INHBA) as well as the glial cell-derived neurotrophic factor (GDNF) and anti-mullerian hormone (AMH). In the data, 5606 unknown transcripts were also expressed (Fig. 2). These matched to genomic sequences of re- lated ruminant species (taurine cattle, yak and goat), but also to domestic pig and ORFV. ORFV is one of the four viruses in the genus Parapoxvirus (PPV) of Poxviridae family, causing contagious ecthyma (orf), primarily in sheep and goat (Crumbie 1998). Recent studies have shown the importance of retrovirus in sheep reproduction and de- velopment (Dunlap et al. 2006). In addition, endogenous retroviruses have been used as genetic markers to ex- plore sheep domestication history (Chessa et al. 2009). Thus, further characterization of those transcripts might provide valuable insights on adaptation of Finnsheep in extreme climate or its exceptional reproductive perfor- mance compared to other breeds. Conclusions Here we present the first transcriptome profiling of the Finnsheep, which is a valuable sheep breed for global sheep farming. In the two studied ovaries from Finnsheep ewes we identified transcripts of 14875 genes. Also, majority of the putative candidate genes related to fertility and out-of-season reproduction were present in both samples. Presence of the majority of putative candidate genes associated with out-of-season breeding in our data not only strengthens previous findings, but also enables further research in this area. Presence of ORFV genes in the data suggests its retro-viral existence in the genome of sheep. Our results showed that ovaries can be a very good candidate tissue for reproductive genomics studies utilizing high throughput sequencing techniques and bioinformatics methods. Authors’ contributions JK, M-HL and JP collected the samples and designed the experiment. KP analyzed the data and wrote the manuscript. JK, M-HL and JP participated in the manuscript writing. All the authors read and approved the final manuscript. Acknowledgements The study was financed by the grants (No. 250633 and 256077; M-HL) from the Academy of Finland and the Na- tional Natural Science Foundation of China (No. 31272413; M-HL). KP’s exchange visit to Prof. Göran Andersson’s lab at SLU, Sweden was supported by the Genomic Resources Research Networking Programme of the European Science Foundation. Eija Korpelainen and Matti Kankainen are acknowledged for their helpful discussions. Com- putational resources were provided by CSC-IT center for Science Ltd. References Anders, S. & Huber, W. 2010. Differential expression analysis for sequence count data. Genome Biology 11: R106. doi:10.1186/ gb-2010-11-10-r106 Anders, S., Pyl, P.T. & Huber, W. 2015. HTSeq - A Python framework to work with high-throughput sequencing data. Bioinformat- ics 31: 166–169. doi:10.1093/bioinformatics/btu638 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 8 Andrews, S. FastQC A Quality Control tool for High Throughput Sequence Data. Retrieved January 18, 2013, from http://www. bioinformatics.babraham.ac.uk/projects/fastqc/ Bonnet, A., Cabau, C., Bouchez, O., Sarry, J., Marsaud, N., Foissac, S., Woloszyn, F., Mulsant, P. & Mandon-Pepin, B. 2013. An over- view of gene expression dynamics during early ovarian folliculogenesis: specificity of follicular compartments and bi-directional dialog. BMC Genomics 14: 904. doi:10.1186/1471-2164-14-904 Chen, L., Liu, K., Zhao, Z., Blair, H.T., Zhang, P., Li, D. & Ma, R.Z. 2012. Identification of sheep ovary genes potentially associated with off-season reproduction. Journal of Genetics and Genomics 39: 181–90. doi:10.1016/j.jgg.2012.03.002 Chessa, B., Pereira, F., Arnaud, F., Amorim, A., Goyache, F., Mainland, I., Kao, Rowland R., Pemberton, J.M., Beraldi, D.,Stear, M.J., Alberti, A.,Pittau, M., Iannuzzi, L., Banabazi, M.H., Kazwala, R.R., Zhang, Y-P., Arranz, J.J., Ali, B.A., Wang, Z., Uzun, M., Dione, M.M., Olsaker, I., Holm, L-E., Saarma, U., Ahmad, S., Marzanov, N., Eythorsdottir, E., Holland, M.J., Ajmone-Marsan, P., Bruford, M.W., Kantanen, J., Spencer, T.E. & Palmarini, M. 2009. Revealing the history of sheep domestication using retrovirus integrations. Sci- ence 324: 532–536. doi:10.1126/science.1170587 Conesa, A., Götz, S., García-Gómez, J.M., Terol, J., Talón, M. & Robles, M. 2005. Blast2GO: a universal tool for annotation, visuali- zation and analysis in functional genomics research. Bioinformatics 21: 3674–3676. doi:10.1093/bioinformatics/bti610 Crumbie, A. 1998. The orf virus: a disease of the farming community. Community Nurse 4: 44–45. Retrieved from http://www. ncbi.nlm.nih.gov/pubmed/9763971 Davis, G.H. 2005. Major genes affecting ovulation rate in sheep. Genetics Selection Evolution 37: 11–24. doi:10.1051/gse Demars, J., Fabre, S., Sarry, J., Rossetti, R., Gilbert, H., Persani, L., Tosser-Klopp, G., Mulsant, P., Nowak, Z., Drobik, W., Martyniuk, E. & Bodin, L. 2013. Genome-wide association studies identify two novel BMP15 mutations responsible for an atypical hyperpro- lificacy phenotype in sheep. PLOS Genetics 9: e1003482. doi:10.1371/journal.pgen.1003482 Dong, J., Albertini, D.F., Nishimori, K., Kumar, T.R., Lu, N. & Matzuk, M.M. 1996. Growth differentiation factor-9 is required during early ovarian folliculogenesis. Nature 383: 531–535. doi:10.1038/383531a0 Dunlap, K.A., Palmarini, M., Varela, M., Burghardt, R.C., Hayashi, K., Farmer, J.L. & Spencer, T.E. 2006. Endogenous retroviruses regulate periimplantation placental growth and differentiation. Proceedings of the National Academy of Sciences of the United States of America 103: 14390–14395. doi:10.1073/pnas.0603836103 Eppig, J. 2001. Oocyte control of ovarian follicular development and function in mammals. Reproduction 122: 829–838. doi:10.1530/ rep.0.1220829 Fahmy, M.H. & Dufour, J.J. 1988. The accumulative effect of Finnsheep breeding in crossbreeding schemes: reproductive perfor- mance. Canadian Journal of Animal Science 68: 69–81. doi:10.4141/cjas88-007 Fogarty, N.M. 2009. A review of the effects of the Booroola gene (FecB) on sheep production. Small Ruminant Research 85: 75– 84. doi:10.1016/j.smallrumres.2009.08.003 Galloway, S.M., McNatty, K.P., Cambridge, L.M., Laitinen, M.P., Juengel, J.L., Jokiranta, T.S., McLaren, R.J., Luiro, K., Dodds, K.G., Montogomery, G.W., Beattie, A.E., Davis, G.H. & Ritvos, O. 2000. Mutations in an oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and infertility in a dosage-sensitive manner. Nature Genetics 25: 279–83. doi:10.1038/77033 Husseneder, C., McGregor, C., Lang, R.P., Collier, R. & Delatte, J. 2012. Transcriptome profiling of female alates and egg-laying queens of the Formosan subterranean termite. Comparative Biochemistry and Physiology. Part D, Genomics & Proteomics 7: 14– 27. doi:10.1016/j.cbd.2011.10.002 Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R. & Salzberg, S. L. 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology 14: R36. doi:10.1186/gb-2013-14-4-r36 Kinsella, R.J., Kähäri, A., Haider, S., Zamora, J., Proctor, G., Spudich, G., Almeida-King, J., Staines, D., Derwent, P., Kerhornou, A., Kersey, P. & Flicek, P. 2011. Ensembl BioMarts: a hub for data retrieval across taxonomic space. Database : The Journal of Biologi- cal Databases and Curation 2011: bar030. doi:10.1093/database/bar030 Koskela, H.L.M., Eldfors, S., Ellonen, P., van Adrichem, A.J., Kuusanmäki, H., Andersson, E.I., Lagström, S., Clemente, M.J., Olson, T., Jalkanen, S.E., Majumder, M.M., Almusa, H., Edgren, H., Lepistö, M., Mattila, P., Guinta, K., Koistinen, P., Kuittinen, T., Penttinen, K., Parsons, A., Knowles, J., Saarela, J., Wennerberg, K., Kallioniemi, O., Porkka, K., Loughran, T.P., Heckman, C.A., Maciejewski, J.P. & Mustjoki, S. 2012. Somatic STAT3 mutations in large granular lymphocytic leukemia. The New England Journal of Medicine 366: 1905–1913. doi:10.1056/NEJMoa1114885 Li, M.-H., Strandén, I., Tiirikka, T., Sevón-Aimonen, M.-L. & Kantanen, J. 2011. A comparison of approaches to estimate the in- breeding coefficient and pairwise relatedness using genomic and pedigree data in a sheep population. PLOS ONE 6: e26256. doi:10.1371/journal.pone.0026256 Maijala, K. 1988. History, recent development and uses of Finnsheep. Journal of Agricultural Science Finland 60: 449–454. Maijala, K. & österberg, S. 1977. Productivity of pure Finnsheep in Finland and abroad. Livestock Production Science 4: 355–377. doi:10.1016/0301-6226(77)90006-9 Mateescu, R.G. & Thonney, M.L. 2010. Genetic mapping of quantitative trait loci for aseasonal reproduction in sheep. Animal Ge- netics 41: 454–9. doi:10.1111/j.1365-2052.2010.02023.x McNatty, K.P., Smith, P., Moore, L.G., Reader, K., Lun, S., Hanrahan, J.P., Groome, N.P., Laitinen, M., Ritvos, O. & Juengel, J.L. 2005. Oocyte-expressed genes affecting ovulation rate. Molecular and Cellular Endocrinology 234: 57–66. doi:10.1016/j.mce.2004.08.013 Mullen, M.P., Hanrahan, J.P., Howard, D.J. & Powell, R. 2013. Investigation of prolific sheep from UK and Ireland for evidence on origin of the mutations in BMP15 (FecX(G), FecX(B)) and GDF9 (FecG(H)) in Belclare and Cambridge sheep. PLOS ONE 8: e53172. doi:10.1371/journal.pone.0053172 Notter, D.R., Cockett, N.E. & Hadfield, T.S. 2003. Evaluation of melatonin receptor 1a as a candidate gene influencing reproduc- tion in an autumn-lambing sheep flock. Journal of Animal Science 81: 912–7. Retrieved from http://www.ncbi.nlm.nih.gov/pub- med/12723079 A G R I C U LT U R A L A N D F O O D S C I E N C E K. Pokharel et al. (2015) 24: 1–9 9 Ramayo-Caldas, Y., Mach, N., Esteve-Codina, A., Corominas, J., Castelló, A., Ballester, M., Estelle, J., Ibanez-Escriche, N., Fernan- dez, A.I., Perez-Enciso, M. & Folch, J.M. 2012. Liver transcriptome profile in pigs with extreme phenotypes of intramuscular fatty acid composition. BMC Genomics 13: 547. doi:10.1186/1471-2164-13-547 Sormunen-Cristian, R. & Suvela, M. 1999. Out-of-season lambing of Finnish Landrace ewes. Small Ruminant Research 31: 265– 272. Retrieved from http://www.sciencedirect.com/science/article/pii/S0921448898001400 Thorvaldsdóttir, H., Robinson, J.T. & Mesirov, J.P. 2013. Integrative Genomics Viewer (IGV): high-performance genomics data visu- alization and exploration. Briefings in Bioinformatics 14: 178–192. doi:10.1093/bib/bbs017 Trapnell, C., Williams, B.A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M.J., Salzberg, S.L., Wold, B.J. & Pachter, L. 2010. Tran- script assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology 28: 511–515. doi:10.1038/nbt.1621 Twine, N.A., Janitz, K., Wilkins, M.R. & Janitz, M. 2011. Whole transcriptome sequencing reveals gene expression and splicing dif- ferences in brain regions affected by Alzheimer’s disease. PLOS ONE 6: e16266. doi:10.1371/journal.pone.0016266 Vinet, A., Drouilhet, L., Bodin, L., Mulsant, P., Fabre, S. & Phocas, F. 2012. Genetic control of multiple births in low ovulating mam- malian species. Mammalian Genome : Official Journal of the International Mammalian Genome Society 23: 727–40. doi:10.1007/ s00335-012-9412-4 Våge, D.I., Husdal, M., Kent, M.P., Klemetsdal, G. & Boman, I.A. 2013. A missense mutation in growth differentiation factor 9 (GDF9) is strongly associated with litter size in sheep. BMC Genetics 14: 1. doi:10.1186/1471-2156-14-1 Xu, Q., Zhao, W., Chen, Y., Tong, Y., Rong, G., Huang, Z., Zhang, Y., Chang, G., Wu, X. & Chen, G. 2013. Transcriptome profiling of the goose (Anser cygnoides) ovaries identify laying and broodiness phenotypes. PLOS ONE 8: e55496. doi:10.1371/journal.pone.0055496