A high resolution atlas of gene expression in the domestic sheep (Ovis aries).

Sheep are a key source of meat, milk and fibre for the global livestock sector, and an important biomedical model. Global analysis of gene expression across multiple tissues has aided genome annotation and supported functional annotation of mammalian genes. We present a large-scale RNA-Seq dataset representing all the major organ systems from adult sheep and from several juvenile, neonatal and prenatal developmental time points. The Ovis aries reference genome (Oar v3.1) includes 27,504 genes (20,921 protein coding), of which 25,350 (19,921 protein coding) had detectable expression in at least one tissue in the sheep gene expression atlas dataset. Network-based cluster analysis of this dataset grouped genes according to their expression pattern. The principle of 'guilt by association' was used to infer the function of uncharacterised genes from their co-expression with genes of known function. We describe the overall transcriptional signatures present in the sheep gene expression atlas and assign those signatures, where possible, to specific cell populations or pathways. The findings are related to innate immunity by focusing on clusters with an immune signature, and to the advantages of cross-breeding by examining the patterns of genes exhibiting the greatest expression differences between purebred and crossbred animals. This high-resolution gene expression atlas for sheep is, to our knowledge, the largest transcriptomic dataset from any livestock species to date. It provides a resource to improve the annotation of the current reference genome for sheep, presenting a model transcriptome for ruminants and insight into gene, cell and tissue function at multiple developmental stages.

Introduction

Sheep (Ovis aries) represent an important livestock species globally and are a key source of animal products including meat, milk and fibre. They remain an essential part of the rural economy in many developed countries and are central to sustainable agriculture in developing countries. They are also an important source of greenhouse gases [1]. Although genetic improvement is often considered to have been less effective in sheep than in the dairy cattle, pig and poultry sectors, advanced genomics-enabled breeding schemes are being implemented in New Zealand and elsewhere [2-4]. A better understanding of functional sequences, including transcribed sequences and the transcriptional control of complex traits such as disease resilience, reproductive capacity and feed conversion efficiency will enable further improvements in productivity with concomitant reductions in environmental impact.

RNA-Sequencing (RNA-Seq) has transformed the analysis of gene expression from the single-gene to the genome-wide scale allowing visualisation of the transcriptome and re-defining how we view the transcriptional control of complex traits (reviewed in [5]). Large-scale gene expression atlas projects have defined the mammalian transcriptome in multiple species, initially using microarrays [6-9] and more recently by the sequencing of full length transcripts or of 5’ ends, for example in the horse [10], and in human and mouse by the FANTOM 5 Consortium [11-13], ENCODE project [14] and Genotype-Tissue Expression (GTEx) Consortium [15].

These efforts have focused mainly on mice and humans, for which there are high quality richly annotated reference genome sequences available as a frame of reference for the identification and analysis of transcribed sequences. Draft reference genome sequences have been established for the major livestock species (chicken, pig, sheep, goat and cattle) over the past decade, yet it is only with the recent deployment of long read sequencing technology that the contiguity of the reference genome sequences for these species has improved. This is exemplified by the recent goat genome assembly [16, 17]. In these species there are still many predicted protein-coding and non-coding genes for which the gene model is incorrect or incomplete, or where there is no informative functional annotation. For example, in the current sheep reference genome, Oar v3.1 (Ensembl release 87) (http://www.ensembl.org/Ovis_aries/Info/Index), 30% of protein-coding genes are identified with an Ensembl placeholder ID [18]. Given the high proportion of such unannotated genes many are likely to be involved in important functions. Large-scale RNA-Seq gene expression datasets can be utilised to annotate and assign function to such unannotated genes [19]. With sufficiently large datasets, genes form co-expression clusters, which can either be generic, associated with a given pathway or be cell-/tissue- specific. This information can then be used to associate a function with genes co-expressed in the same cluster, a logic known as the ‘guilt by association principle’ [20]. Detailed knowledge of the expression pattern can provide a valuable window on likely gene function, as demonstrated in pig [6], sheep [18, 21], human and mouse [8, 9, 22, 23].

A high quality well-annotated reference genome is an exceptionally valuable resource for any livestock species, providing a comparative sequence dataset and a representative set of gene models. The International Sheep Genomics Consortium (ISGC) released a high quality draft sheep genome sequence (Oar v3.1) in 2014 [18]. Included in the sheep genome paper were 83 RNA-Seq libraries from a gestating adult female Texel, 16 day embryo, 8 month old lamb and an adult ram. This Texel RNA-Seq transcriptome significantly improved the annotation of Oar v3.1 and identified numerous genes exhibiting changes in copy number and tissue specific expression [18]. To build on this resource and further improve the functional annotation of Oar v3.1 we have generated a much larger high-resolution transcriptional atlas from a comprehensive set of tissues and cell types from multiple individuals of an outbred cross of two economically important sheep breeds. To maximize heterozygosity we deliberately chose a cross of disparate breeds: the Texel, which is used as a terminal sire as it exhibits enhanced muscling relative to other sheep breeds [24], and the Scottish Blackface, a breed selected for robustness on marginal upland grazing [25].

The sheep gene expression atlas dataset presented here is the largest of its kind from any livestock species to date and includes RNA-Seq libraries from tissues and cells representing all the major organ systems from adult sheep and from several juvenile, neonatal and prenatal developmental time points. Because the tissues were obtained from multiple healthy young adult animals, the atlas may also aid understanding of the function of orthologous human genes. Our aim was to provide a model transcriptome for ruminants and give insight into gene, cell and tissue function and the molecular basis of complex traits. To illustrate the value of the resource, we provide detailed examination of genes implicated in innate immunity and the advantages of cross breeding and provide putative gene names for hundreds of the unannotated genes in Oar v3.1. The entire data set is available in a number of formats to support the research community and will contribute to the Functional Annotation of Animal Genomes (FAANG) project [26, 27].

Results and discussion

Scope of the sheep gene expression atlas dataset

This sheep gene expression atlas dataset expands on the RNA-Seq datasets already available for sheep, merging a new set of 441 RNA-Seq libraries from the Texel x Scottish Blackface (TxBF) with 83 existing libraries from Texel [18]. Details of the new TxBF libraries generated for the sheep gene expression atlas, including the developmental stages sampled, tissue/cell types and sex of the animals are summarised in Table 1. These samples can be grouped into 4 subsets (“Core Atlas”, “GI Tract Time Series”, “Early Development” and “Maternal Reproductive Time Series”). The animals used to generate the four subsets of samples are detailed in S1 Table.

Table 1

Details of the tissues and cell types sequenced to generate the TxBF RNA-Seq dataset for the sheep gene expression atlas.

Subset	Tissue Type	Library Type	Sequencing Depth	Total Number of Libraries	Number of Individuals
Core Atlas	Liver, spleen, ovary, testes, hippocampus, kidney medulla, bicep muscle, reticulum, ileum, thymus, left ventricle	Total RNA-Seq	>100 million paired-end reads per sample	59 (~10 per individual)	3 adult males and 3 adult females
Core Atlas	GI tract, reproductive tract, brain, endocrine, cardiovascular, lymphatic, musculo-skeletal-system, immune cells	mRNA-Seq	>25 million paired-end reads per sample	202 (~45 per individual)	3 adult males and 3 adult females
LPS Time Course(Core Atlas)	BMDM 0h (-LPS) and 7h (+LPS)	Total RNA-Seq	>100 million paired-end reads per sample	12 (2 time points per individual)	3 adult males and 3 adult females
LPS Time Course(Core Atlas)	BMDM 0h, 2h, 4h, 7h, 24h post LPS treatment	mRNA-Seq	>25 million paired-end reads per sample	30 (3 time points per individual)	3 adult males and 3 adult females
GI Tract Time Series	Gastro-intestinal tract	mRNA-Seq	>25 million paired-end reads per sample	108 (~12 tissues per individual)	9 lambs (3 at birth, 3 at one week and 3 at 8 weeks)
Early Development(blastocysts)	Day 7 blastocysts (3 pools of 8)	Nu Gen Single Cell Ovation Kit	>66 million paired-end reads per sample	3	3 pools
Early Development(day 23)	Day 23 whole embryos	Total RNA-Seq	>100 million paired-end reads per sample	3	3 embryos (all male)
Early Development (day 35)	Liver, brain, embryonic fibroblasts	mRNA-Seq	>25 million paired-end reads per sample	7	3 embryos (two female and one male)
Early Development (day 100)	Liver, ovary	mRNA-Seq	>25 million paired-end reads per sample	5	3 female (liver)2 female (ovary)
Maternal Reproductive Time Series (days 23, 35 and 100)	Placenta and ovary	mRNA-Seq	>25 million paired-end reads per sample	12	2 females per time point

Tissues and cells were chosen to cover all major organ systems. All libraries were Illumina 125bp paired end stranded libraries. See S2 Table for a detailed list of the tissues and cell types sequenced.

The “Core Atlas” subset was generated using six adult virgin sheep, approximately 2 years of age. Tissue samples were collected from all major organ systems from 3 males and 3 females to ensure, wherever possible, there were biological replicates from each sex to support an analysis of sex-specific gene expression. In addition, five cell types were sampled, including peripheral blood mononuclear cells (PBMCs) and blood leukocytes. Since macrophages are known to be a highly complex source of novel mRNAs [28], and were not sampled previously, three types of macrophage (+/- stimulation with lipopolysaccharide (LPS)) were included.

For the “GI Tract Time Series” subset of samples we focused on 12 regions of the gastro-intestinal (GI) tract, immediately at birth prior to first feed, at one week and at 8 weeks of age. These time points aimed to capture the transition from milk-feeding to rumination. Embryonic time points were chosen, at days 23, 35 and 100, to detect transcription in the liver, ovary and brain in “Early Development”. Parallel time points were included for placenta and ovary samples from gestating TxBF ewes, comprising the “Maternal Reproductive Time Series” subset. Finally, 3 pools of eight day 7 blastocysts were included to measure transcription pre-implantation and these were also included in the “Early Development” subset.

A detailed list of all tissues and cell types included in each subset of samples can be found in S2 Table. Tissues and cell types were chosen to give as comprehensive a set of organ systems as possible and include those tissues relevant for phenotypic traits such as muscle growth and innate immunity.

Sequencing depth and coverage

Approximately 37x109 sequenced reads were generated from the TxBF libraries, generating approximately 26x109 alignments in total. The raw number of reads and percentage of alignable reads per sample are included in S3 Table. For each tissue a set of expression estimates, as transcripts per million (TPM), were obtained using the high speed transcript quantification tool Kallisto [29]. Kallisto is a new transcriptome-based quantification tool that avoids the considerable bias introduced by the genome alignment step [30]. Gene level expression atlases are available as S1 Dataset and, with expression estimates averaged per tissue per developmental stage, S2 Dataset. We have also made the files containing the expression estimates (S1 Dataset, S2 Dataset and S3 Dataset) available for download through the University of Edinburgh DataShare portal (http://dx.doi.org/10.7488/ds/2112). The data were corrected for library type (as we described in [31] and summarised in S1 Methods). We used Principal Component Analysis (PCA) pre- and post-correction (S1 Fig) for library type to ensure the correction was satisfactory. Hierarchical clustering of the samples is included in (Fig 1) and illustrates both the large diversity and logical clustering of samples included in the dataset.

Fig 1

Hierarchical clustering of the samples included in the sheep gene expression atlas dataset.

Samples of each tissue and cell type from each breed and developmental stage were averaged across individuals for ease of visualisation. The tree was constructed from the Euclidean distances between expression vectors using MEGA v7.0.14 [141] with the neighbour-joining method and edited in the graphical viewer FigTree v1.4.3 [142]. Clustering is biologically meaningful and highlights the lack of any significant effect of library type post-correction. Samples are coloured by organ system.

The O. aries reference genome (Oar v3.1) includes 27,504 loci that are transcribed (20,921 protein coding), of which 25,350 (19,921 protein coding) (95%) were detectable with expression of TPM >1, in at least one tissue from at least one individual, in the sheep gene expression atlas dataset, demonstrating the depth and scope of this dataset. The proportion of genes with detectable expression, after each ‘pass’ with Kallisto (see Methods), is presented in Table 2. After the 'second-pass', only 3% (561) of transcripts from Oar v3.1 (S4 Table) did not meet the minimum detection threshold of TPM > 1 in at least one tissue and therefore were not detected in the sheep atlas dataset. In a minority of cases, genes were missing because they were highly specific to a tissue or cell which was not sampled, such as odontoblasts (which uniquely produce tooth dentin, mediated by DSPP (dentin sialophosphoprotein) [32]). We also did not include any samples taken from the eye which expresses multiple unique proteins e.g., the lens-specific enzyme LGSN (lengsin) [33]. The majority (77%) of the genes not detected in the sheep atlas were unannotated, with no assigned gene name. A small number of these genes (36) lack sequence conservation and coding potential and so are potentially spurious models (S4 Table).

Table 2

The number and percentage of Oar v3.1 protein coding and non-coding genes, with average TPM across all animals > 1 in at least one tissue, in both the TxBF dataset after the Kallisto first and second pass, and after incorporating the existing Texel dataset.

RNA-Seq data used:	Texel x Scottish Blackface Libraries Only					Including the existing Texel dataset
Kallisto index:	First-Pass			Second-Pass (restricted)		Second-Pass (restricted)
Gene type	No. in reference annotation (Oar v3.1)	No. of genes of this type expressed	% genes of this type expressed	No. of genes of this type expressed	% genes of this type expressed	No. of genes of this type expressed	% genes of this type expressed
lincRNA	1858	1548	83.32	0	0	0	0
miRNA	1305	1242	95.17	0	0	0	0
misc RNA	361	310	85.87	0	0	0	0
MT rRNA	2	2	100	0	0	0	0
MT tRNA	22	19	86.36	0	0	0	0
processed pseudogene	43	31	72.09	35	81.40	38	88.37
protein-coding	20921	19921	95.22	20189	96.50	20359	97.31
pseudogene	247	172	69.64	189	76.52	201	81.38
rRNA	305	272	89.18	0	0	0	0
snoRNA	756	717	94.84	0	0	0	0
snRNA	1234	1116	90.44	0	0	0	0
Sum	27054	25350		20413		20598

‘TxBF data’ refers to the present study; ‘Texel data’ is obtained from [18]. The ‘first pass’ Kallisto index contains the known Ovis aries v3.1 cDNAs for both protein-coding and non-protein coding transcripts. The ‘second pass’ Kallisto index is a filtered version of the former, that (a) restricts the RNA space to protein-coding genes, pseudogenes, and processed pseudogenes (so that expression within an equivalent space will be quantified, irrespective of experimental protocol), (b) omits genes that had no detectable expression across all TxBF samples, and (c) includes novel transcript reconstructions further to the de novo assembly of unmapped reads.

Gene annotation

In the Oar v3.1 annotation, 6217 (~30%) of the protein coding genes lack an informative gene name. Whilst the Ensembl annotation will often identify homologues of a sheep gene model, the automated annotation pipeline used is conservative in its assignment of gene names and symbols. Using an annotation pipeline (described in S1 Methods and illustrated in S5 Table) we were able to utilise the sheep gene expression atlas dataset to annotate >1000 of the previously unannotated protein coding genes in Oar v3.1 (S6 Table). These genes were annotated by reference to the NCBI non-redundant (nr) peptide database v77 [34] and assigned a quality category based on reciprocal percentage identity, if any, to one of 9 known ruminant proteomes (S7 Table). A short-list containing a conservative set of gene annotations, to HGNC (HUGO Gene Nomenclature Committee) gene symbols, is included in S8 Table. Many of these genes are found in syntenic regions, and are also supported by the up- and downstream conservation of genes in a related genome, cattle (Bos taurus annotation UMD 3.1). S9 Table contains the full list of genes annotated using this pipeline. Many unannotated genes can be associated with a gene description, but not necessarily an HGNC symbol; these are also listed in S10 Table. We manually validated the assigned gene names on this longer list using network cluster analysis and the “guilt by association” principle.

Network cluster analysis

Network cluster analysis of the sheep gene expression atlas was performed using Miru (Kajeka Ltd, Edinburgh UK), a tool for the visualisation and analysis of network graphs from big data [35-37]. The atlas of unaveraged TPM estimates, available as S1 Dataset, were used for the network cluster analysis. The three blastocyst samples were removed from the network cluster analysis as they were generated using a library preparation method which was not corrected for and created a significant effect of library type. With a Pearson correlation co-efficient threshold of r = 0.75 and MCL (Markov Cluster Algorithm [38]) inflation value of 2.2, the gene-to-gene network comprised 15,129 nodes (transcripts) and 811,213 edges (correlations above the threshold value). This clustering excludes >30% of detected transcripts, most of which had idiosyncratic expression profiles. One of the major sources of unique expression patterns is the use of distinct promoters in different cell types. The transcription factor MITF (Melanogenesis Associated Transcription Factor), for example, does not cluster with any other transcripts in sheep and is known in humans to have at least 7 distinct tissue-specific promoters in different cell types, including macrophages, melanocytes, kidney, heart and retinal pigment epithelium [13].

The resultant correlation network (Fig 2) was very large and highly structured comprising 309 clusters varying in size. Genes found in each cluster are listed in S11 Table and clusters 1 to 50 (numbered in order of size; cluster 1 being made up of 1199 genes) were annotated by hand and assigned a broad functional ‘class’ and ‘sub-class’ (Table 3). Functional classes were assigned based on GO term enrichment [39] for molecular function and biological process (S12 Table) and gene expression pattern, as well as comparison with functional groupings observed in the pig expression atlas [6]. Fig 3 shows a network graph with the nodes collapsed, and the largest clusters numbered 1 to 30, to illustrate the relative number of genes in each cluster and their functional class.

Fig 2

Network visualisation and clustering of the sheep gene expression atlas.

A three-dimensional visualisation of a Pearson correlation gene-to-gene graph of expression levels derived from RNA-Seq data from analysis of sheep tissues and cells. Each node (sphere) in the graph represents a gene and the edges (lines) correspond to correlations between individual measurements above the defined threshold. The graph is comprised of 15,192 nodes (genes) and 811,213 edges (correlations ≥0.75). Co-expressed genes form highly connected complex clusters within the graph. Genes were assigned to groups according to their level of co-expression using the MCL algorithm.

Table 3

Tissue/cell/pathway association of the 50 largest network clusters in the sheep gene expression atlas dataset.

Cluster ID	Number of Transcripts	Profile Description	Class	Sub-Class
1	1199	General	House Keeping	House Keeping (1)
2	987	Brain	Central Nervous System (CNS)	CNS
3	658	Testes–Adult	Reproduction	Gamete Production
4	585	General	House Keeping	House Keeping (2)
5	351	Macrophages	Immune	Macrophages
6	350	Fallopian Tube > Testes	Cilia	Motile Cilia
7	284	Fetal Ovary > Adult Testes	Early Development	Reproduction
8	276	Many Tissues—Highly Variable	Pathway	Cell Cycle
9	265	Fetal Brain > Adult Brain	CNS	CNS
10	247	Skeletal Muscle > Oesophageal Muscle	Musculature	Skeletal Muscle
11	219	Lymph Nodes > Blood > Not Macrophages	Immune	T-Cell and B-Cell
12	215	Thymus > Salivary Gland	Immune	T-Cell
13	186	Fore-Stomachs > Tonsil > Skin	GI Tract	Ruminal Epithelium
14	183	Kidney Cortex > Liver	Renal	Kidney Cortex
15	182	General but not even—highest in muscle	Pathway	Oxidative Phosphorylation
16	158	Epididymis > Vas Deferens	Reproduction	Male
17	153	Liver	Liver	Liver (Hepatocytes)
18	145	Peyer’s Patch, Ileum, Lymph Nodes, Blood	Immune	T-Cell and B-Cell
19	134	Placenta	Gestation	Placental Function
20	119	Epididymis > Testes > Vas Deferens	Reproduction	Male
21	115	General but not even	Pathway	Ribosomal
22	102	Adrenal Gland	Endocrine	Steroid Hormone Biosynthesis
23	102	Placenta	Gestation	Placental Function
24	98	Liver > Small Intestine	Liver	Liver (GI Tract)
25	96	Fetal Liver	Liver	Developing Liver
26	92	Small Intestine > Large Intestine	GI Tract	GI Tract
27	90	Pituitary Gland	Endocrine	Hormone Synthesis
28	85	General highest in reproductive tissues and brain	Cilia	Primary Cilia
29	77	Heart	Musculature	Cardiac Muscle
30	75	Thyroid	Endocrine	Thyroxine Synthesis
31	73	Peyers Patch, Ileum, Lymph Nodes, Blood, Macrophages	Immune	T-Cell and B-Cell
32	69	Salivary Gland, Lymph Node, Blood, Small Intestine	Immune	T-Cell and B-Cell
33	68	Fore-Stomachs Adult—Not Neonates > AMs	GI Tract	Immune
34	65	Small Intestine > Large Intestine	GI Tract	GI Tract
35	61	General—highest in Brain	House Keeping	House Keeping (3)
36	58	Heart Valves	Cardiovascular	Extra Cellular Matrix
37	58	General—highest in Blood	House Keeping	House Keeping (4)
38	56	General—highest in Fetal Brain	House Keeping	House Keeping (5)
39	50	GI Tract—highest in Reticulum	GI Tract	GI Tract
40	49	General—highest in Testes	House Keeping	House Keeping (6)
41	49	General—highest in Ovary	House Keeping	House Keeping (7)
42	47	General—highest in Brain	House Keeping	House Keeping (8)
43	45	Fore-Stomachs > Tonsil > Skin	GI Tract	Ruminal Epithelium
44	45	General—highest in testes	House Keeping	House Keeping (9)
45	44	Macrophage (BMDM + LPS)	Immune	Macrophages (LPS Response TNF)
46	44	Blood, Lymph Nodes	Immune	Blood
47	42	Large Intestine	GI Tract	GI Tract
48	42	General	Pathway	Histones
49	41	Reticulum and Rumen—very variable	GI Tract	Reticulum/Rumen
50	36	Adrenal Gland Medulla	Endocrine	Steroid Hormone Biosynthesis

> indicates decreasing expression profile.

Fig 3

Collapsed node visualisation of the sheep gene expression atlas dataset in two-dimensions to illustrate the relative proportion of genes in each cluster.

Includes 3104 nodes and 138,407 edges with a Pearson correlation value of r = 0.75 and an MCL inflation (MCLi) value of 2.2. Nodes are coloured by tissue/cell type or for broader classes organ system. The largest clusters are numbered from 1 to 30 (see Table 3 for functional annotation). The largest clusters are dominated by either house-keeping genes (1 & 4) or genes associated with transcriptionally rich tissues or cell types, such as brain (2), testes (3) and macrophages (5).

The majority of co-expression clusters included genes exhibiting a specific cell/tissue expression pattern (Fig 4A). There were a few exceptions, including the largest cluster (cluster 1), which contained ubiquitously expressed ‘house-keeping’ genes, encoding proteins that are functional in all cell types. The high proportion of unannotated genes (24% of the 1199 genes) in cluster 1 may reflect the focus of functional genomics on genes exhibiting tissue specific expression, and inferred function in differentiation, leaving those with a house-keeping function uncharacterised [40]. With a few exceptions, the remaining co-expression clusters were composed of genes exhibiting either expression only in a distinct tissue or cell type e.g. macrophages (cluster 5) (Fig 4B (i)) and fetal ovary (cluster 7) (Fig 4B (ii)), or a broader expression pattern associated with a cellular process e.g. oxidative phosphorylation (cluster 15) (Fig 4B (iii)). Some co-expression is shared between two or more organ systems, associated with known shared functions. For example, cluster 15, exhibiting high expression in liver and kidney cortex, is enriched for expression of genes relating to the oxidation-reduction process, transmembrane transport, monocarboxylic acid catabolic process and fatty acid oxidation (S13 Table). It includes numerous genes encoding enzymes involved in amino acid catabolism (e.g. ACY1, AGXT, AGXT2, ASPDH, DPYD, DAO, DDO, EHHADH, HAO1, and HPD) and the rate-limiting enzymes of gluconeogenesis (ALDOB, G6PC, PC, and PCK1). The contributions of kidney and liver to amino acid turnover and gluconeogenesis are well known in humans [41] and rodents [42]. These observations suggest that the shared catabolic pathways of liver and kidney cortex are largely conserved in sheep, but detailed curation of the genes in this cluster could provide further specific insights. Alanine aminotransferase (ALT1; synonym GPT1), which generates alanine from the breakdown of amino acids in muscle and is transported to the liver for gluconeogenesis, is highly expressed in muscle as expected. The glutaminase genes, required for the turnover of glutamine, are absent from tissue or cell type specific clusters; the liver-specific enzyme GLS2 is also expressed in neuronal tissues, as it is in humans [43, 44].

Fig 4

Interrogation of the underlying expression profiles allows regions of the graph to be associated with specific tissues or cell types.

A A three-dimensional visualisation of a Pearson correlation gene to gene network graph (r = 0.75, MCLi = 2.2). Samples of each tissue and cell type from each breed and developmental stage are averaged across individuals for ease of visualisation. Histograms of the averaged expression profile (averaged across individuals for each tissue and cell type for ease of visualisation) of genes in selected clusters are given on the right: B (i) profile of cluster 5 genes whose expression is highest in macrophages; (ii) profile of cluster 7 genes whose expression is highest in fetal ovary and testes; (iii) or a broader expression pattern associated with a cellular process e.g. oxidative phosphorylation (cluster 15). Note that there may be a degree of variation in the expression pattern of individual genes within a cluster which is masked when average profiles are displayed.

The tissue-specific expression patterns observed across clusters showed a high degree of similarity to those observed for pig [6], human and mouse [8, 9]. In some cases we were able to add functional detail to clusters of genes previously observed in pig and human. For example, genes within cluster 6 showed high expression in the fallopian tube and to a lesser extent the testes. Significantly enriched GO terms for cluster 6 included cilium (p = 4.1x10-8), microtubule motor activity (p = 2.9x10-11) and motile cilium (p = 2.8x10-19) suggesting the genes expressed in cluster 6 are likely to have a function related to motile cilia in sperm cells and the fallopian tube. Cluster 6 in the sheep atlas dataset corresponds to cluster 9 in the pig gene expression atlas [6]. Similarly, significantly enriched GO terms for genes in cluster 28 included primary cilium (p = 7.5x10-21) and ciliary basal body (p = 1.9x10-12), indicating the genes in this cluster were associated with the function of primary cilia. Genes within this cluster showed a relatively wide expression pattern, with brain and reproductive tissues exhibiting the highest expression. For both clusters associated with cilial function, significantly enriched GO terms also included cell-cycle related cellular processes and cell-cycle associated genes, supporting the link between cilia and the cell-cycle [45-47].

Cellular processes

The genes within some clusters, rather than being linked to the function of a particular tissue or cell type, showed varying levels of expression across multiple tissues, suggesting their involvement in a universal cellular process (pathway). Significantly enriched GO terms for genes in cluster 8, for example, included ‘cell cycle checkpoint’ (p = 1x10-17) and ‘mitotic cell cycle’ (p = <1x10-30). The variation in expression of these genes across tissues likely reflects variation in the proportion of mitotically active cells. In the same way, it is possible to extract a similar cluster from large cancer gene expression data sets, correlating with their proliferative index [48]. Expression of genes in cluster 15 (n = 182) was detectable in most ovine tissues and cells but with strongly-enriched expression in skeletal and cardiac muscle. The pig gene expression atlas [6] highlighted an oxidative phosphorylation cluster and a mitochondrial/tricarboxylic acid (TCA) cluster. This functional grouping is merged in cluster 15. The majority of the transcripts in cluster 15 are present within the inventory of mitochondrial genes in humans and mice [49], but the reciprocal is not true since many other genes encoding proteins that locate to mitochondria were not found in cluster 15. Many mitochondrial proteins are unrelated to oxidative phosphorylation per se, and are enriched in other tissues including liver and kidney (including mitochondrial enzymes of amino acid and fatty acid catabolism; see above) and intestine.

Cluster 15 also contains several genes associated with myosin and the sarcoplasmic reticulum which may indicate some level of coordination of their function with the oxidative metabolism of glucose. The majority of genes in the corresponding cluster in pig [6] were also present in this cluster with a few notable additions including dihydrolipoamide dehydrogenase (DLD), which encodes a member of the class-I pyridine nucleotide-disulfide oxidoreductase family, and carnitine acetyltransferase (CRAT), a key enzyme in the metabolic pathway in mitochondria [50]. We were able to assign gene names to the following genes associated with oxidative phosphorylation complex I in pig: NDUFA9 (ENSOARG00000009435), NDUFB1 (ENSOARG00000020197), NUDFB8 (assigned to ENSOARG00000015378), and NDUFC2 (ENSOARG00000006694). The gene name PTGES2 (prostaglandin E2 synthase 2) was assigned to ENSOARG00000010878, a gene that was associated with fatty acid (long-chain) beta oxidation in pig [6] but previously unannotated in sheep. Similarly, we assigned the gene name PDHB (pyruvate dehydrogenase E1 component subunit beta), a member of the pyruvate dehydrogenase complex also previously unannotated in sheep, to ENSOARG00000012222. The inclusion of the PDH complex, as well as the mitochondrial pyruvate carriers MPC1 and MPC2, in the muscle-enriched cluster 15 reflects the fact that glucose, giving rise to pyruvate, is the preferred fuel for oxidative metabolism in muscle [51].

By using comparative clustering information in the pig and the “guilt by association” principle we were able to assign with confidence gene names and putative function to the majority of unannotated genes in clusters 8 (cell-cycle) and 15 (oxidative phosphorylation) (S13 Table). Expression of cell cycle and metabolic genes has recently been shown to be positively correlated with dry matter intake, ruminal short chain fatty acid concentrations and methane production in sheep [52]. In the same study a weak correlation between lipid/oxo-acid metabolism genes and methane yield was also identified suggesting that the unannotated genes in these clusters are likely to be relevant in addressing methane production in ruminants [52].

The GI tract

Stringent coexpression clustering requires that each transcript is quantified in a sufficiently large number of different states to establish a strong correlation with all other transcripts with which it shares coordinated transcription and, by implication, a shared function or pathway. The impact of this approach was evident from the pig gene expression atlas [6] which was effective at dissecting region-specific gene expression in the GI tract. We have generated a comparable dataset for sheep. In ruminants, the rumen, reticulum and omasum are covered exclusively with stratified squamous epithelium similar to that observed in the tonsil [18, 21]. Each of these organs has a very distinctive mucosal structure, which varies according to region sampled [53]. A network cluster analysis of regions of the GI tract from sheep has been published [21] using the Texel RNA-Seq dataset [18]. These co-expression clusters are better refined in this larger atlas, because many of the genes that are region-specific in the GI tract are also expressed elsewhere. We have, in addition, expanded the available dataset for the GI tract to include samples from neonatal and juvenile lambs.

The postnatal development of the sheep GI tract is of particular interest because of the pre-ruminant to ruminant transition, which occurs over 8 weeks from birth. Genes in cluster 33 showed low levels of expression in neonatal lambs and a gradual increase into adulthood. Enriched GO terms for this cluster include regulation of interleukin 6 (IL6) production (p = 0.0016) and keratinocyte differentiation (p = 1.7x10-8) (S12 Table). The cluster includes genes such as HMGCS2, HMGCL and BDH1, required for ketogenesis, an essential function of rumen metabolism, as well as CA1 (carbonic anhydrase 1), implicated in the rumen-specific uptake of short chain fatty acids. The cluster does not contain any of the solute carriers implicated in nutrient uptake in the rumen, suggesting that these are more widely-expressed and/or regulated from multiple promoters [21]. The only carrier that is rumen-enriched is SLC9A3 (also known as NHE3), the key Na-H antiporter previously implicated in rumen sodium transport in both sheep and cattle [54]. Other genes in cluster 33, for example, IL36A and IL36B, are thought to influence skin inflammatory response by acting directly on keratinocytes and macrophages and indirectly on T-lymphocytes to drive tissue infiltration, cell maturation and cell proliferation [55]. Many of these genes might also be part of the acute phase immune response, by regulating production of key cytokines such as IL-6 and thus mediating activation of the NF-κB signaling pathways. Nuclear factor (NF)-κB and inhibitor of NF-κB kinase (IKK) proteins regulate innate- and adaptive-immune responses and inflammation (reviewed in [56]). Expression of many of these genes is likely to change as the immune system develops which we will describe in detail in a dedicated network cluster analysis of the GI tract developmental time series dataset. The genes in this cluster therefore appear to be involved both in the onset of rumination and in innate immunity (which could be associated with the population of the rumen microbiome).

Innate and acquired immunity

Several clusters exhibited a strong immune signature. Clusters 11, 12, 18, 31 and 32, for example, contained genes with a strong T-lymphocyte signature [57] with high levels of expression in immune cell types and lymphoid tissues. Significantly enriched GO terms for cluster 12, for example, included T-cell differentiation (p = 2.10x10-11), immune response (p = 0.00182) and regulation of T-cell activation (p = 0.00019) (S12 Table). Manual gene annotation using Ensembl IDs within this cluster revealed the majority were T-cell receptors and T-cell surface glycoproteins (S14 Table). Interestingly, ENSOARG00000008993 represents a gene with no orthologues to other species within the Ensembl database, but partial blast hits to T-lymphocyte surface antigen Ly-9 in mouflon, goat, bison and buffalo in the NCBI database. The ‘true’ gene LY9, a member of the signalling lymphocyte activation molecule (SLAM) family [58], is also unannotated in sheep and is assigned to ENSOARG00000008981, having multiple orthologues in other placental mammals. We have assigned ENSOARG00000008993 the gene name ‘LY9-like’ and the symbol LY9L, and suggest this gene plays a role in T-lymphocyte pathogen recognition.

Other immune clusters exhibited a macrophage-specific signature, with subsets highly expressed in alveolar macrophages (AMs), monocyte derived macrophages (MDMs) and bone marrow derived macrophages (BMDMs) (cluster 5) and two defined clusters of genes induced in BMDMs stimulated with LPS (cluster 45 and 52). Known macrophage-specific surface markers, receptors and proinflammatory cytokines predominated in these clusters, in addition to numerous unannotated genes, with as yet undefined but probable immune function (S15 Table). For example, the CD63 antigen, which mediates signal transduction events, was assigned to ENSOARG00000011313 and BST2 (bone marrow stromal cell antigen 2), which is part of the interferon (IFN) alpha/beta signaling pathway, to ENSOARG00000016787. A third cluster of LPS-inducible genes in macrophages, cluster 64, contained a subset of the IFN-inducible antiviral effector genes, including DDX58, IFIT1, IFIT2, MX1, MX2, RSAD2 and XAF1, which are induced in mouse and humans through the MyD88-independent TLR4 signaling pathway via autocrine IFNB1 signaling (reviewed in [59]). Many other components of this pathway identified in LPS-stimulated human macrophages [60] were either not annotated, or not clustered, and will be the target of detailed annotation efforts in the macrophage dataset.

Significantly enriched GO terms for the macrophage-specific cluster 5 included ‘response to lipopolysaccharide’ (p = 7.2x10-7), and ‘toll-like receptor signaling pathway’ (p = 3.2x10-5). Many of the genes in this cluster are known components of the innate immune response in mammals. Interleukin 27 (IL-27), is a heterodimeric cytokine which has pro- and anti-inflammatory properties and a diverse effect on immune cells including the regulation of T-helper cell development, stimulation of cytotoxic T-cell activity and suppression of T-cell proliferation [61]. ADGRE1 encodes the protein EGF-like module-containing mucin-like hormone receptor-like 1 (EMR1; also known as F4/80), a classic macrophage marker in mice [62]. Several genes in cluster 5 encode proteins exclusively expressed in macrophages and monocytes. One such gene, CD163, encodes a member of the scavenger receptor cysteine-rich (SRCR) superfamily, which protects against oxidative damage by the clearance and endocytosis of hemoglobin/haptoglobin complexes by macrophages, and may also function as an innate immune sensor of bacteria [63].

One of the largest macrophage populations in the body occupies the lamina propria of the small and large intestine [64]. They are so numerous that the expression of macrophage-related genes can be detected within the total mRNA from intestine samples. As noted previously in the pig, one can infer from the expression profiles that certain genes that are highly-expressed in AMs are repressed in the intestinal wall [6]. We proposed that such genes, which included many c-type lectins and other receptors involved in bacterial recognition, were necessary for the elimination of inhaled pathogens, where such responses would be undesirable in the gut [6]. In the sheep, there was no large cohort of receptors that showed elevated expression in AMs relative to MDMs or BMDMs, and that were absent in the gut wall. Only a small cluster (115) of 13 genes showed that profile, including the phagocytic receptor VSIG4 (CRiG), which is a known strong negative regulator of T-cell proliferation and IL2 production [65] and SCIMP, recently identified as a novel adaptor of Toll-like receptor signaling that amplifies inflammatory cytokine production [66]. Six previously unannotated genes within this small cluster included the E3 ubiquitin ligase, MARCH1, and likely members of the paired immunoglobulin type receptor and SIGLEC families, which cannot be definitively assigned as orthologues.

Interestingly, macrophage colony-stimulating factor receptor (CSF1R), which controls the survival, proliferation and differentiation of macrophage lineage cells [67, 68], was not within the macrophage-specific cluster 5. Instead, CSF1R was in a small cluster (cluster 102) along with several other macrophage-specific genes including the C1Q complex. As in humans and mice [9, 12, 69], CSF1R was also expressed in sheep placenta. In humans and mice, placental (trophoblast) expression is directed from a separate promoter [11]. The small number of genes co-expressed with CSF1R are likely either co-expressed by trophoblasts as well as macrophages (as is C1Q in humans; see BioGPS (http://biogps.org/dataset/GSE1133/geneatlas-u133a-gcrma/) [9, 70]), or highly-expressed in placenta-associated macrophages.

Early development and reproduction

The sheep gene expression atlas dataset includes multiple libraries from early developmental time points. Three of the larger clusters of co-expressed genes showed high levels of expression in the fetal ovary (cluster 7), fetal brain (cluster 9) and fetal liver (cluster 25). ‘Testis-specific’ genes, particularly those involved in meiosis and gametogenesis, might also be expressed in the fetal ovary undergoing oogenesis [71, 72]. Our dataset from sheep appears to validate this hypothesis, since genes within cluster 7 exhibited higher levels of expression in the fetal ovary and to a lesser extent the testes. Several genes were expressed both in the testes and the fetal ovary, including testis and ovary specific PAZ domain containing 1 (TOPAZ1), which has been shown in sheep to be expressed in adult male testes and in females during fetal development with a peak during prophase I of meiosis [73], and fetal and adult testis expressed 1 (FATE1), which is strongly expressed in spermatogonia, primary spermatocytes, and Sertoli cells in seminiferous tubules in mouse and humans [74]. Significantly enriched GO terms for genes within cluster 7 included ‘female gonad development’ (p = 4.9x10-6), ‘spermatogenesis’ (p = 4.6x10-8) and ‘growth factor activity’ (p = 5x10-5) (S12 Table).

Several important genes for embryonic development were also co-expressed in cluster 7. The germ-cell specific gene SRY-box 30 (SOX30) encodes a member of the SOX (SRY-related HMG-box) family of transcription factors involved in the regulation of embryonic development and in the determination of cell fate [75]. Growth differentiation factor 3 (GDF3) encodes a protein required for normal ocular and skeletal development. Although it is a major stem cell marker gene [76], it has not previously been linked to germ cell expression. Similarly, POU class 5 homeobox 1 (POU5F1) encodes a transcription factor containing a POU homeodomain that controls embryonic development and stem cell pluripotency [76] but is also required for primordial germ cell survival [77]. The expression of these genes in tissues containing germ cells in sheep suggests they contribute to meiosis and cellular differentiation. These observations illustrate the utility of the sheep as a non-human model for the study of gametogenesis.

Cluster 7 also includes two related oocyte-derived members of the transforming growth factor-β (TGFB1) superfamily, growth differentiation factor 9 (GDF9) and bone morphogenetic protein 15 (BMP15), which are essential for ovarian follicular growth and have been shown to regulate ovulation rate and influence fecundity in sheep [78, 79]. Lambing rate is an important production trait in sheep and can vary between breeds based on single nucleotide polymorphism (SNP) mutations in key genes influencing ovulation rate (reviewed in [79, 80]). A number of the known fertility genes in sheep (reviewed in [81, 82]), such as the estrogen receptor (ESR) and the Lacaune gene (B4GALNT2) were not present in cluster 7, which may be because they are not expressed in the ovary at the time points chosen for this study. Detailed analysis of the expression of key genes during early development in the fetal ovary in comparison with the ovary from the adult and gestating ewes may provide additional insights.

Sex specific differences in gene expression

Sex-specific differences in gene expression have been reported in humans [83, 84] mice [85, 86], cattle [87, 88] and pigs [89, 90]. We examined male and female biased gene expression in the sheep atlas dataset by calculating the average TPM per sex for each gene and the female:male expression ratio (S3 Dataset). Twenty genes exhibited strongly sex biased expression (S16 Table); 13 were female-enriched and 7 were male-enriched. Among the male enriched genes was thyroid stimulating hormone beta (TSHB), which is expressed in thyrotroph cells in the pituitary gland and part of a neuro-endocrine signaling cascade in sheep [91]. Expression of TSHB in the pituitary gland of male TxBF was 3.6-fold higher than in female TxBF sheep. A similar sex bias has been observed in rats in which males exhibit significantly higher TSHB expression in the pituitary gland than females [92].

Other genes exhibiting similarly large sex specific fold-changes included keratin 36 (KRT36) which was expressed 6.6-fold higher in the reticulum of male relative to female sheep and VSIG1 (V-Set and immunoglobulin domain containing 1), which is required for the differentiation of glandular gastric epithelia [93]. VSIG1 showed 4-fold greater expression in the female pylorus relative to the male. The unannotated gene ENSOARG00000020792 exhibited large fold change in male biased expression in immune tissues including popliteal and prescapular lymph node, tonsil and Peyer’s patch. This gene has a detectable blast hit to “immunoglobulin kappa-4 light chain variable region” and is a 1:1 orthologue with an unannotated gene in cattle, ENSBTAG00000045514, with ≥70% reciprocal percentage identity and conservation of synteny. The dN/dS for ENSOARG00000020792 suggests it is evolving rapidly (dN/dS > 2). Male-biased genes are known to evolve quickly, as are immune genes [94]. GO term enrichment for the set of genes with five-fold sex-biased expression in at least one TxBF tissue (S17 Table) revealed that the genes enriched in females were predominately involved with the immune response while genes enriched in the male were broadly associated with muscle and connective tissue. This is likely to reflect inherent differences between the two sexes in allocation of resources towards growth or immunity. Genes exhibiting sex specific expression might therefore be relevant in sexual diamorphism in disease susceptibility, for example. Additionally, pregnancy is characterized by significant and complex changes in immune parameters which is likely to impact on sex specific differences in gene expression. Further investigation is warranted into the complexity of sex specific differences in gene expression throughout development.

Differential expression of genes between the Texel and TxBF

The majority of commercially-produced livestock are a cross of two or more different production breeds with distinct desired traits [95]. For example, in the UK, the crossing of lighter upland sheep breeds with heavier lowland meat breeds optimises carcass quality, lambing rate, growth rate and survivability [95]. In developing countries, sustainable crossing of indigenous small ruminants with elite western breeds is one approach to improve productivity [96, 97]. An RNA-Seq dataset of this size from an outbred cross of two disparate sheep breeds provides an opportunity to investigate differential gene expression in a purebred parental line and crossbred animals. We compared gene expression across tissues in the F1 crossbred (TxBF) animals (generated by crossing Texel rams with Scottish Blackface ewes; Fig 5A) with the purebred Texel animals included in the previous sheep gene expression atlas dataset [18].

Fig 5

Upregulation of genes in the crossbred TxBF relative to the purebred Texel.

A A Texel sire was crossed with Scottish Blackface dam to create the F1 Texel x Scottish Blackface individuals. B The top 20 genes showing the greatest up-regulation (as absolute fold change) between the crossbred TxBF and purebred Texel individuals. Genes associated with the function of skeletal muscle and connective tissue are indicated.

A gene was considered differentially expressed (DE) between the purebred Texel and hybrid TxBF if (a) it was expressed at ≥1 TPM in both Texel and TxBF (considering TPM to be the mean of all replicates per tissue), (b) the fold change (ratio of TxBF TPM to Texel TPM) was ≥2 in ≥25% of the tissues in which expression was detected (stipulating no minimum number of tissues, but noting that 23 tissues are common to Texel and TxBF), and (c) the fold change was ≥5 in at least 1 tissue. Fold changes of all genes expressed at ≥1 TPM in both breeds are given in S18 Table. The GO terms enriched in the set of DE genes (n = 772) with higher expression in the TxBF than the Texel were predominantly related either to muscle or brain function (S19 Table). The top 20 genes showing the largest up-regulation (shown as absolute fold-change) in the TxBF relative to the purebred Texel sheep are illustrated in Fig 5B. Enriched molecular function GO terms for the set of genes differentially expressed between TxBF and Texel sheep include ‘iron ion binding’ (p = 3.6x10-4), and ‘cytoskeletal protein binding’ (p = 7.8x10-6), biological process terms include ‘cellular iron ion homeostasis’ (p = 6.3x10-4) and cellular component terms include ‘sarcomere’ (p = 6.6x10-7) and ‘collagen trimer’ (p = 5.5x10-7) (S19 Table).

Numerous genes with structural, motor and regulatory functions were highly expressed in TxBF compared to Texel bicep muscle, with approximately 5- to 18-fold expression increases for various members of the collagen (COL1A1, COL1A2, COL3A1) and myosin families (MYH2, MYL4), along with CSRP3 (a mechanosensor) [98], FMOD (fibromodulin, a regulator of fibrillogenesis) [99], keratocan (KERA, a proteoglycan involved in myoblast differentiation) [100], matrix metalloproteinase 2 (MMP2, a proteolytic enzyme associated with muscle regeneration) [101], and calsequestrin 2 (CASQ2, one of the most abundant Ca2+-binding proteins in the sarcoplasmic reticulum, essential for muscle contraction) [102].

Genes enriched in muscle are of particular biological and commercial interest because Texel sheep exhibit enhanced muscling and less fat [103], due to a single nucleotide polymorphism (SNP) in the 3’ untranslated region of the myostatin gene MSTN (synonym GDF-8) which generates an illegitimate miRNA binding site resulting in translational inhibition of myostatin synthesis and contributing to muscular hypertrophy [24, 104]. Because heterozygotes have an intermediate phenotype, cross breeding of homozygous mutant Texel sheep with animals homozygous for the normal allele transmits something of the Texel muscle phenotype to the offspring. An effect on muscle synthesis in the TxBF animals could be related to the myostatin genotype; genes with higher expression in the Texel than in the cross may be targets for myostatin inhibition [105], while those with lower expression in the Texel than in the cross may be directly or indirectly activated by myostatin and hence involved in the cessation of muscle differentiation. In cattle the myostatin mutation is associated with the downregulation of collagen genes including COL1A1 and COL1A2 [105]. This is consistent with the observation that these genes have higher expression in the heterozygous TxBF animals than the Texel animals. Since myostatin also regulates muscle fibre type [106] by suppressing the formation of fast-twitch fibres, individuals homozygous for inactivating myostatin mutations are likely to exhibit increased fast-twitch fibres [107]. Many of the genes up-regulated in the TxBF relative to the Texel animals (e.g. CSRP3 and CASQ2) are known to be specifically expressed in slow-twitch muscle [106, 108], and several down-regulated genes are associated with fast-twitch muscle (e.g. TNNC2, TNNI2 and SERCA1) [109, 110]. Consequently, the difference between the cross-breed and pure Texel is in part attributable to an increased contribution of slow-twitch fibres, which in turn has been associated with desirable meat quality traits [111] highlighting the potential advantages of cross-breeding.

Enriched GO terms related to brain function include the ‘myelin sheath’ (p = 6.1x10-8) and the ‘internode region of the axon’ (p = 5.2x10-5) (S19 Table). Candidate genes of particular interest were expressed in the cerebellum (S18 Table). For instance, in the TxBF relative to the Texel animals, there were approximately 8-fold expression increases in cochlin (COCH, which regulates intraocular pressure) [112] and brevican (BCAN, which functions throughout brain development in both cell-cell and cell-matrix interactions) [113, 114] and a 10-fold expression increase for myelin-associated oligodendrocyte basic protein, MOBP, which was previously unannotated in sheep (ENSOARG00000002491) and has a function in late-stage myelin sheath compaction [115]. A 15-fold expression increase was observed for oligodendrocytic paranodal loop protein (OPALIN, a transmembrane protein unique to the myelin sheath [116]) and a 10-fold increase for another unannotated gene myelin basic protein (MBP), which we have assigned to ENSOARG00000004374 [117]. Although these examples of a neuroendocrine-specific effect of cross-breeding are speculative, they are of interest as Scottish Blackface sheep exhibit both improved neonatal behavioural development [118] and more extensive foraging behaviour than lowland breeds such as the Texel, travelling further distances, covering greater areas and exploring higher altitudes [119].

Visualisation of the expression atlas

We have provided the TxBF sheep gene expression atlas as a searchable database in the gene annotation portal BioGPS (http://biogps.org/dataset/BDS_00015/sheep-atlas/). By searching the dataset via the following link (http://biogps.org/sheepatlas/) the expression profile of any given gene can be viewed across tissues. An example profile of the myostatin (MSTN) gene from sheep is included in Fig 6. BioGPS allows comparison of expression profiles across species and links to gene information for each gene [70, 120, 121]. The Sheep Atlas BioGPS expression profiles are based on TPM estimates from the alignment-free Kallisto output for the TxBF libraries, averaged across samples from each developmental stage for ease of visualization (S2 Dataset). It is important to note that there may be a degree of variation in the expression pattern of specific genes between individuals which is masked when the average profiles are displayed. In addition, to allow comparison between species BioGPS requires each gene have an Entrez ID, which is not the case for all genes in Oar v3.1 and as a consequence these genes do not have visualisable expression profiles in BioGPS. The expression profiles of the genes without Entrez IDs can be found in S1 Dataset and S2 Dataset.

Fig 6

Screenshot of the expression profile of the sheep myostatin (MSTN) gene within the BioGPS online portal.

Expression estimates from the TxBF sheep gene expression atlas dataset are available via the BioGPS database (http://biogps.org/dataset/BDS_00015/sheep-atlas/). This provides a searchable database of genes, with expression profiles across tissues and cells for each gene displayed as histograms via the following link, http://biogps.org/sheepatlas/. The BioGPS platform supports searching for genes with a similar profile, allows access to the raw data, and links to external resources. It also provides the potential for comparative analysis across species, for example with the expression profiles for pig.

In parallel to the alignment-free Kallisto method, we also used an alignment-based approach to RNA-Seq processing, with the HISAT aligner [122] and StringTie assembler [123] (detailed in S1 Methods). These alignments will be published as tracks in the Ensembl genome browser in the short term and integrated into the next Ensembl genome release for sheep.

Conclusions

This work describes the transcriptional landscape of the sheep across all major organs and multiple developmental stages, providing the largest gene expression dataset from a livestock species to date. The diversity of samples included in the sheep transcriptional atlas is the greatest from any mammalian species, including humans. Livestock provide an attractive alternative to rodents as models for human conditions in that they are more human-like in their size, physiology and transcriptional regulation, as well as being economically important in their own right. Non-human models are required to study the fundamental biology of healthy adult mammals and as such this dataset represents a considerable new resource for understanding the biology of mammalian tissues and cells.

In this sheep transcriptional atlas gene expression was quantified at the gene level across a comprehensive set of tissues and cell-types, providing a starting point for assigning function based on cellular expression patterns. We have provided functional annotation for hundreds of genes that previously had no meaningful gene name using co-expression patterns across tissues and cells. Future analysis of this dataset will use the co-expression clusters to link gene expression to observable phenotypes by highlighting the expression patterns of candidate genes associated with specific traits from classical genetic linkage studies or genome-wide association studies (GWAS). Gene expression datasets have been used in this way to characterise cell populations in mouse [22, 23] and in the biological discovery of candidate genes for key traits in sheep [124-126] and pigs [127, 128]. We have already utilised the dataset to examine the expression patterns of a set of candidate genes linked to mastitis resistance [129] in sheep, including comparative analysis with a recently available RNA-Seq dataset from sheep lactating mammary gland and milk samples [130]. The research community will now be able to use the sheep gene expression atlas dataset to examine the expression patterns of their genes or systems of interest, to answer many of the outstanding questions in ruminant biology, health, welfare and production.

Improving the functional annotation of livestock genomes is critical for biological discovery and in linking genotype to phenotype. The Functional Annotation of Animal Genomes Consortium (FAANG) aims to establish a data-sharing and research infrastructure capable of efficiently analysing genome wide functional data for animal species [26, 27]. This analysis is undertaken on a large scale, including partner institutions from across the globe, to further our understanding of how variation in gene sequences and functional components shapes phenotypic diversity. Analysis of these data will improve our understanding of the link between genotype and phenotype, contribute to biological discovery of genes underlying complex traits and allow the development and exploitation of improved models for predicting complex phenotypes from sequence information. The sheep expression atlas is a major asset to genome assembly and functional annotation and provides the framework for interpretation of the relationship between genotype and phenotype in ruminants.

Methods

Animals

Approval was obtained from The Roslin Institute’s and the University of Edinburgh’s Protocols and Ethics Committees. All animal work was carried out under the regulations of the Animals (Scientific Procedures) Act 1986. Three male and three female Texel x Scottish Blackface sheep of approximately two years of age were acquired locally and housed indoors for a 7–10 day “settling-in period” prior to being euthanased (electrocution and exsanguination). It was not recorded whether the females were in estrus. Nine Texel x Scottish Blackface lambs were born at Dryden Farm Large Animal Unit, produced by mating Texel rams (4 different sires were used in total) with Scottish Blackface ewes. Three neonatal lambs were observed at parturition and euthanised immediately prior to their first feed, three lambs were euthanised at one week of age prior to rumination (no grass was present in their GI tract) and three at 8 weeks of age once rumination was established. The lambs were euthanised by schedule one cranial bolt. To obtain developmental tissues six Scottish Blackface ewes were each mated to a Texel ram (one of the four different sires used to produce the lambs) and scanned to ensure successful pregnancy at 21 days. Two were euthanised at 23 days, two at 35 days and two at 100 days gestation (electrocution followed by exsanguination). Corresponding time points (day 23, 35 and 100) from gestating TxBF ewes mated with a Texel ram (one of four different sires as above), euthanised with the same method as the Blackface ewes, were also included. All the animals were fed ad libitum on a diet of hay and 16% sheep concentrate nut, with the exception of the lambs pre-weaning who suckled milk from their mothers. Details of the animals sampled are included in S1 Table.

Tissue collection

Tissues (95 tissues/female and 93 tissues/male) and 5 cell types were collected from three male and three female adult Texel x Scottish Blackface (TxBF) sheep at two years of age. The same tissues were collected from nine lambs, 3 at birth, 3 at one week and 3 at 8 weeks of age. Three embryonic time points were also included: three day-23 TxBF whole embryos, three TxBF day 35 embryos from which tissue was collected from each region of the basic body plan and three day 100 TxBF embryos from which 80 tissues were collected. Reproductive tissue from the corresponding time points from 6 TxBF ewes, 2 at each gestational time point, was also collected. In addition, 3 pools of 8 day 7 blastocysts from abattoir derived oocytes (of unknown breed) fertilized with Texel semen were created using IVF.

The majority of tissue samples were collected into RNAlater (AM7021; Thermo Fisher Scientific, Waltham, USA) and a subset were snap frozen, including lipid rich tissues such as adipose and brain. To maintain RNA integrity all tissue samples were harvested within an hour from the time of death. A detailed list of the tissues collected and sequenced can be found in S2 Table. Within the scope of the project we could not generate sequence data from all the samples collected and have archived the remainder for future analysis. Sample metadata, conforming to the FAANG Consortium metadata standards, for all the samples collected for the sheep gene expression atlas project has been deposited in the BioSamples database under project identifier GSB-718 (https://www.ebi.ac.uk/biosamples/groups/SAMEG317052).

Isolation of cell types

All cell types were isolated on the day of euthanasia. Bone marrow cells were isolated from 10 posterior ribs as detailed for pig [131]. BMDMs were obtained by culturing bone marrow cells for 7 days in complete medium: RPMI 1640, Glutamax supplement (35050–61; Invitrogen, Paisley, UK), 20% sheep serum (S2263; Sigma Aldrich, Dorset, UK), penicillin/streptomycin (15140; Invitrogen) and in the presence of recombinant human CSF-1 (rhCSF-1: 104 U/ml; a gift of Chiron, Emeryville, CA) on 100-mm2 sterile petri dishes, essentially as described previously for pig [131]. For LPS stimulation the resulting macrophages were detached by vigorous washing with medium using a syringe and 18-g needle, washed, counted, and seeded in tissue culture plates at 106 cells/ml in CSF-1–containing medium. The cells were treated with LPS from Salmonella enterica serotype minnesota Re 595 (L9764; Sigma-Aldrich) at a final concentration of 100 ng/ml as previously described in pig [131] and harvested into TRIzol® (15596018; Thermo Fisher Scientific) at 0, 2, 4, 7 and 24 h post LPS treatment before storing at -80°C for downstream RNA extraction.

PBMCs were isolated as described for pig [132]. MDMs were obtained by culturing PBMCs for 7 days in CSF-1–containing medium, as described above for BMDMs, and harvesting into TRIzol® (15596018; Thermo Fisher Scientific). Alveolar macrophages were obtained by broncho-alveolar lavage of the excised lungs with 500ml sterile PBS (Mg2+ Ca2+ free) (P5493; Sigma Aldrich). The cells were kept on ice until processing. To remove surfactant and debris cells were filtered through 100uM cell strainers and centrifuged at 400 × g for 10 min. The supernatant was removed and 5ml red blood cell lysis buffer (420301; BioLegend, San Diego, USA) added to the pellet for 5 min; then the cells were washed in PBS (Mg2+ Ca2+ free) (P5493; Sigma Aldrich) and centrifuged at 400 × g for 10 min. The pellet was collected, resuspended in sterile PBS (Mg2+ Ca2+ free) (P5493; Sigma Aldrich), and counted. Alveolar macrophages were seeded in 6-well tissue culture plates in 2ml complete medium: RPMI 1640, Glutamax supplement (35050–61; Invitrogen), 20% sheep serum (S2263; Sigma Aldrich), penicillin/streptomycin (15140; Invitrogen) in the presence of rhCSF1 (104 U/ml) overnight.

Blood leukocytes were isolated as described in [133]. Whole blood was spun at 500 x g for 10 min (no brake) to separate the buffy coats. These were then lysed in ammonium chloride lysis buffer (150mM NH4Cl, 10mM NaHCO3, 0.1mM EDTA) for 10 min on a shaking platform, then centrifuged at 4°C for 5 min at 500 x g. The resultant blood leukocyte pellets were stored in 1ml of RNAlater (AM7021; Thermo Fisher Scientific) at -80°C.

To isolate embryonic fibroblasts we harvested a day 35 embryo whole and transferred to outgrowth media: DMEM, high glucose, glutamine, pyruvate (11995065; Thermo Fisher Scientific), FBS (Fetal Bovine Serum) (10500056; Thermo Fisher Scientifc), MEM NEAA (11140068; Thermo Fisher Scientific), penicillin/streptomycin (15140; Invitrogen), Fungizone (15290018; Amphotericin B; Thermo Fisher Scientific), Gentamicin (15750037; Thermo Fisher Scientific). In a sterile flow hood the head was removed and the body cavity eviscerated. The remaining tissue was washed 3 times in PBS (Mg2+ Ca2+ free) (P5493; Sigma Aldrich) with penicillin/streptomycin (15140; Invitrogen). 5ml of Trypsin-EDTA solution (T4049; Sigma Aldrich) was added and the sample incubated at 37°C for 5 min then vortexed and incubated for an additional 5min at 37°C. 3ml of solution was removed and filtered through a 100uM cell strainer, 5ml of outgrowth media was then passed through the strainer and combined with the sample. The sample was centrifuged at 200 x g for 3 min and the pellet resuspended in 9ml of outgrowth media before splitting the sample between 3x T75 flasks. The process was then repeated for the remaining 2ml of sample left after the digestion with Trypsin (T4049; Sigma Aldrich). Embryonic fibroblasts were incubated for 5–7 days (until 80–90% confluent) then harvested into TRIzol® (15596018; Thermo Fisher Scientific).

RNA extraction and library preparation

RNA was extracted using the same method as the Roslin RNA-Seq samples included in the sheep genome project detailed in [18]. For each RNA extraction <100mg of tissue was processed. Care was taken to ensure snap frozen samples remained frozen prior to homogenisation, and any cutting down to the appropriate size was carried out over dry ice. Tissue samples were first homogenised in 1ml of TRIzol (15596018; Thermo Fisher Scientific) with either CK14 (432–3751; VWR, Radnor, USA) or CKMIX (431–0170; VWR) tissue homogenising ceramic beads on a Precellys Tissue Homogeniser (Bertin Instruments; Montigny-le-Bretonneux, France). Homogenisation conditions were optimised for tissue type but most frequently 5000 rpm for 20 sec. Cell samples which had previously been collected in TRIzol (15596018; Thermo Fisher Scientific) were mixed by pipetting to homogenise. Homogenised (cell/tissue) samples were then incubated at room temperature for 5 min to allow complete dissociation of the nucleoprotein complex, 200μl BCP (1-bromo-3-chloropropane) (B9673; Sigma Aldrich) was added, then the sample was shaken vigorously for 15 sec and incubated at room temperature for 3 min. The sample was centrifuged for 15 min at 12,000 x g, at 4°C to separate the homogentate into a clear upper aqueous layer (containing RNA), an interphase and red lower organic layers (containing the DNA and proteins), for three min. DNA and trace phenol was removed using the RNeasy Mini Kit (74106; Qiagen Hilden, Germany) column purification, following the manufacturer’s instructions (RNeasy Mini Kit Protocol: Purification of Total RNA from Animal Tissues, from step 5 onwards). RNA quantity was measured using a Qubit RNA BR Assay kit (Q10210; Thermo Fisher Scientific) and RNA integrity estimated on an Agilent 2200 Tapestation System (Agilent Genomics, Santa Clara, USA) using the RNA Screentape (5067–5576; Agilent Genomics) to ensure RNA quality was of RINe > 7.

RNA-Seq libraries were prepared by Edinburgh Genomics (Edinburgh Genomics, Edinburgh, UK) and run on the Illumina HiSeq 2500 sequencing platform (Illumina, San Diego, USA). Details of the libraries generated can be found in S2 Table. Libraries were sequenced at a depth of either >100 million, >60 million or >25 million paired end reads per sample depending upon to which subset of samples they belonged. In each case depth refers to the number of paired end reads, therefore a depth of >100 million reads represents ~100M forward + 100M reverse. A subset of 10 tissue samples and BMDMs at 0 h and 7h (+/-LPS) (Table 1), from each individual, were sequenced at a depth of >100 million strand-specific 125bp paired-end reads per sample using the standard Illumina TruSeq total RNA library preparation protocol (Ilumina; Part: 15031048, Revision E). These samples were chosen to include the majority of transcriptional output, as in [134]. An additional 40 samples from the tissues and cell types collected per individual (44/female and 42/male), were selected and sequenced at a depth of >25 million strand-specific paired-end reads per sample using the standard Illumina TruSeq mRNA library preparation protocol (poly-A selected) (Ilumina; Part: 15031047 Revision E).

In addition to the samples from the 6 adults, tissue was also collected from other developmental time points. The GI tract tissues collected from the 9 TxBF lambs, 3 at birth, 3 at one week of age and 3 at 8 weeks of age were sequenced at a depth of >25 million paired-end reads per sample using the Illumina mRNA TruSeq library preparation protocol (poly-A selected) as above. Of the early developmental time points, the three 23 day old embryos from TxBF sheep were sequenced at >100 million paired-end reads using the Illumina total RNA TruSeq library preparation protocol (as above), while the other embryonic samples and the ovary and placenta from the gestating ewes were sequenced at a depth of >25 million paired-end reads per sample using the Illumina mRNA TruSeq library preparation protocol (as above). In addition, three libraries were generated using the NuGen Ovation Single Cell RNA-Seq System (0342-32-NUG; NuGen, San Carlos, USA) from pooled samples of 8 blastocysts (as in [135]), and sequenced at a depth of >60 million paired-end reads per sample. A detailed list of prepared libraries, including library type can be found in S2 Table.

To identify spurious samples we used sample-to-sample correlation, of the transposed data from S1 Dataset, in Miru (Kajeka Ltd, Edinburgh, UK) [37]. The sample-to-sample graph is presented in S2 Fig. The expression profiles of any samples clustering unexpectedly (i.e. those not found within clusters of samples of the same type/biological replicate) were examined in detail. Generally the correlation between samples was high, although 4 spurious samples, 4 sets of swapped samples, and 3 samples where a sample collection issue (collection of muscle rather than oesophageal tissue) had occurred were identified. These samples were either relabeled or removed as appropriate and are listed in S20 Table.

Data quality control and processing

The raw data, in the form of.fastq files, for the 438 TxBF libraries is deposited in the European Nucleotide Archive under study accession number PRJEB19199 (http://www.ebi.ac.uk/ena/data/view/PRJEB19199)). A description of the abbreviations used for each sequencing library in submission PRJEB19199 is included in S21 Table. The data submission to the ENA includes experimental metadata prepared according to the FAANG Consortium metadata and data sharing standards. The RNA-Seq data processing methodology and pipelines are described in detail in S1 Methods. For each tissue a set of expression estimates, as transcripts per million (TPM), were obtained using the high speed transcript quantification tool Kallisto [29]. In total, the expression atlas utilised approximately 26 billion (pseudo) alignments (S3 Table), capturing a large proportion of protein-coding genes per tissue (S22 Table).

The accuracy of Kallisto is dependent on a high quality index (reference transcriptome) [29], so in order to ensure an accurate set of gene expression estimates we employed a ‘two-pass’ approach. We first ran Kallisto on all samples using as its index the Oar v3.1 reference transcriptome. We then parsed the resulting data to revise this index. This was in order to include, in the second index, those transcripts that should have been there but were not (i.e. where the reference annotation is incomplete), and to remove those transcripts that should not be there but were (i.e. where the reference annotation is poor quality and a spurious model has been introduced). For the first criterion, we obtained the subset of reads that Kallisto could not align, assembled those de novo into putative transcripts (S1 Methods), then retained each transcript only if it could be robustly annotated (by, for instance, encoding a protein similar to one of known function) and showed coding potential (S23 Table). For the second criterion, we identified those members of the reference transcriptome for which no evidence of expression could be found in any of the hundreds of samples comprising the atlas. These were then discarded from the index. Finally, this revised index was used for a second iteration of Kallisto, generating higher-confidence expression level estimates. This improved the capture rate of protein-coding genes (S24 Table). A detailed description of this process can be found in S1 Methods.

We complemented this alignment-free method with a conventional alignment-based approach to RNA-Seq processing, using the HISAT aligner [122] and StringTie assembler [123]. A detailed description of this pipeline is included in S1 Methods. This assembly is highly accurate with respect to the existing (Oar v3.1) annotation, precisely reconstructing almost all exon (96%), transcript (98%) and gene (99%) models (S25 Table). Although this validates the set of transcripts used to generate the Kallisto index, we did not use HISAT/StringTie to quantify expression. This is because a standardised RNA space is necessary to compare data from mRNA-Seq and total RNA-Seq libraries [31], which cannot be applied if expression is quantified via genomic alignment. Unlike alignment-free methods, however, HISAT/StringTie can be used to identify novel transcript models (S26 Table), particularly for ncRNAs, which will be described in detail in a dedicated analysis. We will publish the alignments from HISAT and StringTie as tracks in the Ensembl genome browser in the short term and integrate the alignments into Ensembl and Biomart in the next Ensembl release for sheep.

Inclusion of additional RNA-Seq datasets from Sheep

Additional RNA-Seq data was obtained from a previous characterisation of the transcriptome of 3 Texel sheep included in the release of the current sheep genome Oar v3.1 [18]. The dataset included tissues from an adult Texel ram (n = 29), an adult Texel ewe (n = 25) and their female (8–9 month old) lamb (n = 28), plus a whole embryo (day 15 gestation) from the same ram-ewe pairing. The raw read data from the 83 Texel samples incorporated into this dataset and previously published in [18] is located in the European Nucleotide Archive (ENA) study accession PRJEB6169 (http://www.ebi.ac.uk/ena/data/view/PRJEB6169)). The metadata for these individuals is included in the BioSamples database under Project Identifier GSB-1451 (https://www.ebi.ac.uk/biosamples/groups/SAMEG317052)). A small proportion of the tissues included in the Texel RNA-Seq dataset were not sampled in the TxBF gene expression atlas. Those unique to the Texel are largely drawn from the female reproductive, integument and nervous systems: cervix, corpus luteum, ovarian follicles, hypothalamus, brain stem, omentum and skin (side and back). Details of the Texel RNA-Seq libraries including tissue and cell type are included in S27 Table. The Texel samples were all prepared using the Illumina TruSeq stranded total RNA protocol with the Ribo-Zero Gold option for both cytoplasmic and mitochondrial rRNA removal, and sequenced using the Illumina HiSeq 2500 (151bp paired-end reads) [18]. As above, Kallisto was used to estimate expression level for all samples, using the revised reference transcriptome (from the ‘second pass’) as its index.

Correcting for the effect of multiple library types

To correct for the confounding effect of multiple library types we applied a batch effect correction. We have previously validated this method using a subset of the sheep expression atlas samples from BMDMs (+/- LPS) sequenced both as mRNA and total RNA libraries [31]. As described above, for the Kallisto second pass, we constrained the Kallisto index to contain only the transcripts of protein-coding genes, pseudogenes and processed pseudogenes, the majority of which are poly(A)+ and so are present in both mRNA-Seq and total RNA-Seq samples. We then calculated, per gene, the ratio of mean TPM across all mRNA-Seq libraries to mean TPM across all total RNA-Seq libraries. Given the scope of the tissues sampled for both library types (all major organ systems from both sexes and from different developmental states), neither mean is likely to be skewed by any tissue-specificity of expression. As such, any deviations of this ratio from 1 will reflect variance introduced by library type/depth. Thus, to correct each gene’s set of expression estimates for this effect of library type, we multiplied all total RNA-Seq TPMs by this ratio. To validate this approach we used principal component bi-plot analysis, described and shown in S1 Methods and S1 Fig.

Gene expression, network cluster analysis and annotation

Network cluster analysis of the sheep gene expression atlas was performed using Miru (Kajeka Ltd, Edinburgh, UK) [35-37]. In brief, similarities between individual gene expression profiles were determined by calculating a Pearson correlation matrix for both gene-to-gene and sample-to-sample comparisons, and filtering to remove relationships where r < 0.75. A network graph was constructed by connecting the remaining nodes (genes) with edges (where the correlation exceeded the threshold value). This graph was interpreted by applying the Markov Cluster algorithm (MCL) [38] at an inflation value (which determines cluster granularity) of 2.2. The local structure of the graph was then examined visually. Genes with robust co-expression patterns, implying related functions, clustered together, forming cliques of highly interconnected nodes. A principle of ‘guilt by association’ was then applied, i.e. the function of an unannotated gene could be inferred from the genes it clustered with [20, 136]. Expression profiles for each cluster were examined in detail to understand the significance of each cluster in the context of the biology of sheep tissues and cells. Clusters 1 to 50 (Table 2) were assigned a functional ‘class’ and ‘sub-class’ manually by first determining if multiple genes within a cluster shared a similar biological function based on both gene ontology [39], determined using the Bioconductor package ‘topGO’ [137] (GO term enrichment for clusters 1 to 50 is shown in S12 Table). We then compared the clusters with tissue- and cell-specific clusters in other large-scale network-based gene expression analyses including the pig gene expression atlas [6], the human protein atlas [69, 72, 138] and the mouse atlas [9, 139, 140]. More specific annotation of the GI tract clusters in sheep was based on network and pathway analysis from the sheep genome paper and a subsequent satellite publication [18, 21]. The gene component of all clusters can be found in S11 Table.

We assigned gene names to unnannotated genes in Oar v3.1 based on their co-expression pattern, tissue specificity, and reciprocal percent identity to a set of nine known ruminant proteomes (S7 Table). The annotation pipeline is described in detail in S1 Methods and included a set of quality categories summarised in S5 Table. We were able to assign gene names to >1000 previously unannotated genes in Oar v3.1. Candidate gene names are given as both a shortlist (S8 Table) and a longlist (S9 Table), the latter intended for manual review as informative annotations may still be made without every one of the above criteria being met.

Principal component analysis with all samples plotted in two dimensions using their projections onto the first two principal components.

For each sample, each gene’s expression level is taken as the mean TPM across all replicates, before (A) and after (B) any batch effect correction. Samples are coloured by organ system. Ellipses indicate confidence intervals of 95%. The shape of each point indicates each sample’s library type: mRNA-seq (circle) or total RNA-seq (triangle). Blastocyst samples are excluded for clarity as they are generated using a different experimental protocol. Before correction, points can be partitioned by shape (to the left and right of sub-Fig A), suggesting a batch effect–variation introduced by library type confounds variation by tissue type. After correction (sub-Fig B), there is no notable axis of variation that partitions points by shape–consequently, variation introduced by library type (a batch effect) does not confound variation by tissue type (which is biologically meaningful).

(TIFF)

Click here for additional data file.

Sample-to-sample network graph analysis was used to validate the samples included in the sheep gene expression atlas dataset.

Each node represents a sample and each edge its connectivity to other samples in the dataset. A correlation of r = 0.75 split the graph into 10 different clusters. The largest cluster (cluster 1) included the majority of samples (‘mixed tissues’), most of which were transcriptionally similar, while the remainder of the clusters comprised samples with distinctive transcriptional signatures such as macrophages (3) and abomasum. Spurious samples were easily identified if they were present in a cluster comprised of samples from a different tissue or cell type. Pearson correlation r = 0.75, MCLi = 2.2, nodes = 481 and edges = 23,903.

(TIFF)

Click here for additional data file.

Additional methods.

(DOCX)

Click here for additional data file.

Gene expression level atlas as transcripts per million (unaveraged).

(TXT)

Click here for additional data file.

Gene expression level atlas as transcripts per million (averaged across biological replicates for each developmental stage).

(TXT)

Click here for additional data file.

Sex-biased expression atlas (based on expression estimates from 3 adult male and 3 adult female TxBF sheep).

(TXT)

Click here for additional data file.

Overview of animals used to generate the sheep atlas tissue subsets.

(XLSX)

Click here for additional data file.

Details of library type and tissue/cell samples used in each subset of samples.

(XLSX)

Click here for additional data file.

Number of reads, and number of aligned reads, per sample.

(XLSX)

Click here for additional data file.

Genes undetected (TPM < 1) in every tissue/cell line of the expression atlas.

(XLSX)

Click here for additional data file.

Quality categories for automated gene annotations.

(XLSX)

Click here for additional data file.

Proportion of unannotated protein-coding genes assigned probable gene names.

(XLSX)

Click here for additional data file.

Source of ruminant proteome data.

(XLSX)

Click here for additional data file.

Candidate gene names for unannotated Oar v3.1 protein-coding genes: shortlist.

(XLSX)

Click here for additional data file.

Candidate gene names for unannotated Oar v3.1 protein-coding genes: longlist.

(XLSX)

Click here for additional data file.

Candidate gene descriptions for all unannotated Oar v3.1. protein-coding genes (potentially informative in the absence of a gene name).

(XLSX)

Click here for additional data file.

Genes within each co-expression cluster.

(XLSX)

Click here for additional data file.

GO term enrichment for co-expression clusters 1 to 50.

(XLSX)

Click here for additional data file.

Manual annotation of cluster 15 (genes involved in oxidative phosphorylation).

(XLSX)

Click here for additional data file.

Manual annotation of unannotated genes in cluster 12 (genes with a T-cell signature).

(XLSX)

Click here for additional data file.

Manual annotation of unannotated genes in cluster 5 (genes with an alveolar macrophage signature).

(XLSX)

Click here for additional data file.

Genes with five-fold sex-biased expression in at least one TxBF tissue.

(XLSX)

Click here for additional data file.

GO term enrichment for the set of genes with five-fold sex-biased expression in at least one TxBF tissue.

(XLSX)

Click here for additional data file.

Fold changes in expression level between TxBF and Texel sheep.

(XLSX)

Click here for additional data file.

GO term enrichment for the set of genes differentially expressed between TxBF and Texel sheep.

(XLSX)

Click here for additional data file.

Details of the removed, re-labelled and swapped TxBF libraries.

(XLSX)

Click here for additional data file.

Details of the abbreviations used for each sequencing library submitted to the European Nucleotide Archive under study accession number PRJEB19199.

(XLSX)

Click here for additional data file.

Number of genes with detectable expression, per tissue.

(XLSX)

Click here for additional data file.

Coding potential of putative novel CDS.

(XLSX)

Click here for additional data file.

Number of genes with detectable expression, per gene type.

(XLSX)

Click here for additional data file.

Proportion of Oar v3.1 gene, exon and transcript models in the StringTie assembly.

(XLSX)

Click here for additional data file.

Novel transcript models in the StringTie assembly.

(XLSX)

Click here for additional data file.

Details of the additional RNA-Seq libraries included from Texel sheep [18].

(XLSX)

Click here for additional data file.