Effect of Gender and Residual Feed Intake Phenotype on the Transcriptome of Liver and Skeletal Muscle Tissue in Purebred Simmental Beef Cattle

Background: The selection of cattle with enhanced feed eciency is of paramount importance with regard to reducing feed costs in the beef industry. However the role, if any, of gender to the underlying molecular control of feed eciency in cattle is not currently known. Global transcriptome proling was undertaken on liver and skeletal muscle biopsies from Simmental heifers and bulls divergent in residual feed intake (RFI) feed eciency phenotype, in order to identify differentially expressed genes that may be associated with this trait. Results: We identied 5 genes (p<0.001; false discovery rate (FDR) <0.1) to be differentially expressed in skeletal muscle between high and low RFI heifers with all 5 transcripts being up-regulated in the low RFI phenotype. Among these differentially expressed genes, all transcripts were involved in oxidative phosphorylation and mitochondrial homeostasis. A total of 11 genes (p<0.001;FDR <0.1) were differentially expressed in hepatic tissue between high and low RFI bulls with 8 transcripts being up-regulated and 3 being down-regulated in the low RFI phenotype. These differentially expressed genes were related to oxidative response, protein mediation and cell signalling. No genes were identied as differentially expressed in both heifer liver and bull muscle analyses. Conclusions: Results from this study show a clear effect of gender to the underlying molecular control of RFI in cattle, which may be attributable to differences in the physiological age between heifers and bulls. Despite this we have highlighted a number of genes that may hold potential as molecular biomarkers for RFI cattle.


Background
Global agriculture is currently faced with the ambitious challenge of feeding a rapidly increasing global population, expected to peak at 9.2 billion by 2075 [1]. This necessary increase in agricultural outputs must also be achieved within the current con nes of arable land availability, thus it is essential that animal production systems become more e cient for the continued sustainability of the beef production sector. Animal feed can account for up to 75% of the variable costs in beef production systems; hence reduction of these costs is of paramount importance [2]. Feed e ciency (dietary nutrient utilisation) in beef cattle is a trait of major economic importance [3]. Indeed, there is signi cant phenotypic and genetic variation among beef animals in their ability to convert dietary derived nutrients into saleable product [4,5,6]. Thus, by improving feed e ciency it is possible to reduce feed intake in cattle while still maintaining growth and skeletal muscle gain and ultimately contributing to beef production pro tability and sustainability.
Residual feed intake (RFI), rst described by Koch et al. (1963) [7] is one such measure of feed e ciency, which can be de ned as the difference between an animal's actual versus its predicted feed intake based on average daily gain (ADG) and metabolic weight [8]. Through measuring an animal's inherent RFI, feede cient animals which have a low-RFI value and consume less than expected as well as feed ine cient animals which have a high-RFI and consume more than expected value may be identi ed [3,[9][10][11][12][13][14][15].
Indeed, it has been shown that growing beef cattle divergent for RFI can consume up to 20% less feed than their counterparts for the same level of performance [3,10,11,[13][14][15][16][17][18][19]. This coupled with a moderate heritability estimate for RFI [20] provides a feasible method for effective utilisation of this trait in production systems through genomic selection processes. However, although moderately heritable, a challenge remains to reliably and cost-effectively identify feed-e cient, low-RFI animals and to proliferate their genetics through animals breeding programmes. For example, the primary impediment to genetic progress and adoption of selection strategies based on RFI is both the large-scale logistics and expense of measuring individual animal feed intake and body weight.
In order to overcome the aforementioned limitations, studies have sought to uncover the underlying biology governing the trait, with the twin-goal for the identi cation of molecular biomarkers [21]. This not only provides an attractive alternative to direct measurement of dietary intake on large numbers of animals [6], but it also allows for a better understanding of the biological mechanisms underlying RFI which is paramount for progressing genomic selection. However although results in the literature have reported roles for biological processes including lipogenesis and the immune system toward variation in RFI primarily through transcriptional pro ling, there is a distinct lack of commonality of key genes contributing to RFI across these studies. This is undoubtedly due to the multi-faceted nature of RFI as well as the contribution of factors including breed, gender, stage of development as well as production management and dietary intake test period length to the trait [21]. For example a recent study that evaluated the molecular control of cattle divergent in RFI across three contrasting breed types, reported only 5 genes as commonly differentially expressed between high and low-RFI groups [22]. Furthermore, of these 5 genes only one, SCD, was differentially expressed in only one other molecular based RFI evaluation [23]. This lack of commonality, although due to confounding experimental designs represents a major short-fall toward the progress of genomics selection. Equally there is a dearth of information in relation to the potential effect of gender on subsequent variation in RFI, particular under similar rearing conditions. Currently it is not known whether the same molecular mechanisms and genes are contributing to variation in RFI in heifers and bulls destined for beef production. Thus the objective of this study was to evaluate the underlying biology regulating RFI through transcriptome sequencing in Simmental heifers and bulls reared under the same conditions from birth. It is estimated that two thirds of the variation in RFI is due to variation in resting energy expenditure [5], with both muscle and liver representing important metabolic tissues with skeletal muscle accounting for over 50% body weight [12] and liver accounting for 18-25% of total oxygen body consumption [24]. Thus our efforts were focused towards examining the transcriptional alterations of these tissues between beef heifers and bulls divergent for high and low-RFI.

Animal performance
Details of live weight gain, feed intake, and animal performance are presented in Tables 1 and 2 for bulls and heifers, respectively. Brie y, bulls and heifers had a mean ADG of 1.8 kg and 1.3 kg and dry matter intake (DMI) of 9.3 kg and 9.2 kg respectively during the RFI measurement period. Residual feed intake averaged − 0.04 and 0.04 and ranged from − 1.92 to 1.04 kg of dry matter(DM)/d and − 0.43 to 0.68 of DM/d for bulls and heifers respectively. Bulls and heifers ranked as high RFI consumed 10 and 15% more than their low RFI counterparts (p < 0.05), respectively. Bulls of high and low RFI did not differ (p > 0.05) in initial bodyweight (BW), nal BW or ADG and likewise for heifers.  3 Back fat change is mean of fat depth at end of intake trial -mean of fat depth at start of intake trial. approximately 85% of reads aligned to the bovine reference genome (UMD3.1) across both tissues analysed. A total of 5 genes (p < 0.001; false discovery rate (FDR) < 0.1) were identi ed to be differentially expressed in skeletal muscle tissue between high and low RFI heifers with all 5 transcripts being upregulated in the low RFI phenotype. Figure 1 provides a visual representation of the correlation matrix pertaining to the differentially expressed genes (DEGs) from heifer skeletal muscle tissue, indicating all DEGs to be highly correlated with each other. A total of 11 genes (p < 0.001; FDR < 0.1) were identi ed to be differentially expressed in hepatic tissue between high and low RFI bulls with 8 transcripts being upregulated and 3 being down-regulated in the low RFI phenotype. Figure 2 provides a visual representation of the correlation matrix pertaining to the DEGs from bull hepatic tissue, highlighting the most correlated genes. Further details of these DEGs are provided in Tables 3 and 4. Correlation analyses between DEGs and RFI for each heifer and bull group, identi ed mostly signi cant associations. For example DEGs from the heifer skeletal muscle analysis were all signi cantly negatively associated (p < 0.05) with RFI, with the exception of cytb, which tended towards signi cance (p = 0.0622). Similarly within the liver analyses of bulls divergent for RFI, 7 of the 11 DEGs were signi cantly associated with RFI (p < 0.05), with DBP, MANF and GMPPB tending towards signi cance (p < 0.1). ACTA2 expression was not signi cantly associated with RFI in the liver tissue of bulls (p = 0.3078). Correlation results between DEGs and RFI value for each group are presented in Table 5. Among the (approximately) 7,000 genes detected in heifer hepatic tissue, none were found to be differentially expressed. Similarly, no differential expression was detected between RFI phenotypes in the (approximately) 6,000 genes expressed in skeletal muscle tissue from bulls.

Differentially Expressed Genes And Pathway Analysis
Functional analysis using both the up-and down-regulated genes was performed using DAVID software and revealed that in hepatic tissue from bulls the DEG were related to metabolism of amino sugars and nucleotides, posttranslational modi cation, protein turnover and chaperones and protein processing in the endoplasmic reticulum (Table 6). Functional analysis using DAVID in heifer skeletal muscle indicated the DE genes to be involved in oxidative phosphorylation and mitochondrial homeostasis ( Table 7). All DEG were successfully mapped to a molecular or biological pathway and/or category in the Ingenuity Pathway analysis (IPA) database. DEG were analysed and separated according to their biological function within IPA. The top canonical pathways affected by RFI in heifer skeletal muscle tissue and hepatic tissue of bulls are presented in Table 8. These included enriched pathways related to mitochondrial function and oxidative phosphorylation in heifers and Aldosterone Signalling, GDPmannose Biosynthesis, Nuclear factor erythroid 2-related factor 2 (NRF2) mediated oxidative stress response, and Eukaryotic Initiation Factor 2 (EIF2) signalling in bulls.

Discussion
Residual feed intake is a prime feed-e ciency trait to target within beef production systems; this is not only due to the moderate heritability of the trait but also due to its independence from production traits used to calculate it [21]. Furthermore, data from our own group has also shown that RFI is a repeatable trait in beef cattle [11]. However although well suited toward genetic propagation through genomic selection, results on the underlying genes governing the trait are con icting across studies [21,25]. This is undoubtedly due to the multifaceted nature of the trait as well as the in uence of both animal and management factors toward the RFI phenotype. Additional individual experimental parameters may also be contributing to the lack of consistency across results, including for example the measurement of the RFI phenotype and the length of the dietary intake test period [26]. Furthermore RFI calculated within a speci c population, as is the case in RFI transcriptional pro ling studies, only re ects natural variation within that speci c population, thus the level of divergence may be quite variable across different studies. In addition to these limitations, genomic selection models do not differentiation on the basis of gender, thus if molecular biomarkers are to be successfully employed for a trait it is essential to determine whether biological processes and speci c key genes are regulating economically important traits such as RFI across differing gender types. Thus again the objective of this study was to evaluate any potential effects of differing gender to the underlying biological mechanisms regulating variation in RFI in both Simmental heifers and bulls. The animals used in the current study were from a purebred, wellcharacterised herd, reared as a contemporary group from birth with similar genetics, thus permitting a more equitable comparison of the effect of phenotypic RFI ranking. Our evaluations were focused towards both the liver and skeletal muscle tissue, the metabolic activities of which are both essential for overall body homeostasis and e ciency of an animal with both organs being highly abundant in mitochondria [27]. Muscle accounts for approximately 50% of body mass and 25% of basal metabolic activity of an animal and plays an important role in resting energy expenditure [28]. Additionally the liver is a highly oxidative organ accounting for 18-25% of total oxygen body consumption that is responsible for metabolising lipids, proteins and carbohydrates into biologically useful molecules [24]. We hypothesised that due to the metabolic importance of these organs, that variation in feed e ciency and energy expenditure (measured here using RFI) is likely to be re ected in the transcriptome of tissue from these organs. Understanding the essential biological processes contributing to variation in RFI is critical to elucidating the genetic basis, for this trait.
Although the RFI values pertaining to the animals used in this study showed clear signi cant divergence between high and low-RFI groups in both heifers and bulls, we failed to identify any DEGs within the liver tissue of the heifers and the skeletal muscle tissue of the bulls. Similarly, other studies evaluating the molecular control of RFI divergence in crossbred steers reported no DEGs following correction for multiple testing [29][30][31]. Therefore the main observation of the current study is the inconsistent effect of RFI across (i) gender and (ii) tissue for DEG pro les despite a 10% and 15% difference in DMI between high and low RFI heifers and bulls, respectively, with no difference in ADG in the current study. While both genders were of similar age and were reared under the same conditions, it is apparent that both genders were at distinct physiological stages at the time of sample collection. Similarly in a targeted gene expression study, complementary to this current study, we observed a signi cant effect of gender in the expression of lipogenesis genes within the subcutaneous adipose tissue, this is despite there being no signi cant difference in subcutaneous fat measurements [32]. Typically heifers display earlier physiological maturity when compared to bulls of the same age, thus the identi cation of a gender effect in an energy storage accretion tissue such as adipose tissue was perhaps unsurprising, however the potential effect of physiological age and stage of maturity on metabolic tissues at the same age in the current study was potentially unexpected. Our results indicate that it will be unlikely that the key genes derived from the tissues examined in this study will be accurate predictors of genetic potential for RFI across gender. However this does not mean that other tissues or organs within the body may provide more comparable results and thus reliable biomarkers, based on animal age and not physiological age or stage of development. Despite this, we determined possible molecular mechanisms and biological functions in uencing RFI in beef cattle within tissue and gender subgroups. These data highlight a relationship between RFI and the transcriptomic networks involved in mitochondrial function in the skeletal muscle of heifers and evidence for a relationship between RFI and the aldosterone signalling pathway and the NRF2 mediated oxidative stress pathway in hepatic tissue of bulls.
We identi ed 5 genes that were statistically signi cantly differentially expressed in skeletal muscle of high and low RFI heifers. Although we did not identify large numbers of DEGs, the level of transcriptional differences are consistent with the ndings of other previously published RFI based studies [33,34,35]. Among the 5 DEGs, all 5, COX1 (cytochrome c oxidase subunit 1), ND5 (NADH-dehydrogenase 5), ND6 (NADH-dehydrogenase 6), CYTB (cytochrome b), and COX3 (cytochrome c oxidase subunit 3), were upregulated in low RFI heifers compared to their high RFI counterparts. Additionally all of these genes were signi cantly associated with RFI, with the exception of cytb, which only tended towards a signi cant association, further indicated a role for these genes to RFI variation in the skeletal muscle tissue of heifers. Interestingly all of these genes are components of the electron transport chain (ETC) in the mitochondrial inner membrane and are of major importance to overall energetic e ciency. ND5 and ND6 are subunits of the enzyme complex NADH dehydrogenase (ubiquinone) or complex I of the ETC. CYTB is the main subunit of coenzyme Q: cytochrome c-oxidoreductase or complex III of the ETC. COX1 is the main subunit of cytochrome c oxidase or otherwise known as complex IV of the ETC and COX3 is a transmembrane subunit of this same complex. Previous transcriptomic experiments have reported a relationship between DEG related to the complexes of ETC and RFI [12,36]. For example, increased levels of COXII (complex IV) and NADH dehydrogenase subunits have been shown to be associated with e cient animals at the protein and transcriptome level [12,[36][37][38][39] in more e cient animals. However, although genes of the ETC have previously been reported in relation to variation in RFI, this is the rst report of an up-regulation of COX1, ND5, ND6, CYTB and COX3 in skeletal muscle tissue of cattle of low vs. high RFI. The lack of commonality of speci c key genes underlying processes such as ETC to RFI phenotype may be due to the differences in breed types employed across studies, with Simmental cattle utilised in the current study. For example, in a study examining the effect of various breeds on the molecular control of RFI in liver tissue, Mukiibi et al. [22] recorded only 5 genes as commonly differentially expressed across three differing breed types, however an evaluation of the biological processes showed a clear commonality for the underlying biological control of RFI irrespective of breed, but with different key genes dependent on the breed [22].
Using gene ontology (GO) enrichment analysis through the DAVID platform, we identi ed important processes underlying feed e ciency variation in skeletal muscle of heifers. The most enriched cluster contained a plethora of functions related to mitochondrial metabolism, in particular oxidative phosphorylation. Similarly, the top canonical pathways identi ed by IPA were mitochondrial function and oxidative phosphorylation. Taken together, these analyses highlight a greater capacity for mitochondrial function in the low RFI animals. The relationship between mitochondria and RFI has been addressed previously and mitochondrial dysfunction and oxidative stress have been implicated as contributing to variation in feed e ciency across varying species [12,36,[40][41][42][43]. Mitochondria are highly dynamic organelles that are responsible for 90% of the energy production in the body and are major reactive oxygen species (ROS) regulators [44]. It seems likely that variation in mitochondrial function could contribute to variation in energy utilisation. Furthermore, it has been demonstrated that feed e cient animals exhibit greater capacity to modulate conditions of oxidative stress [45].
Feed e cient animals have been shown to have a higher activity of all enzymes of the ETC across multiple species including broilers and lambs [40,41,43] and the present study focused on cattle is in agreement with this. The DEG in the current study encode proteins involved in Complexes I, III and IV of the ETC indicating an impaired oxidative phosphorylation system in the skeletal muscle of the less e cient heifers. These results complement previous research suggesting an association between decreased respiration capacity and increased ROS production in less e cient animals [41,46]. Moreover, the results of the present study are reinforced by the observation of Kong et al. [39] in which the mitoproteome was skewed towards high feed e ciency birds despite no difference in mitochondrial DNA between phenotypes, suggesting an increase in mitochondrial activity in the high feed e ciency phenotype [39].
Within the hepatic tissue of bulls we identi ed 11 DEGs, 8 of which were up-regulated in the low RFI animals. These included; HSPA5 (78 kDa glucose -regulated protein precursor), CRELD2 (cysteine rich with EGF like domains 2), HYOU1 (hypoxia up-regulated protein 1 precursor), MANF (mesencephalic astrocyte derived neurotrophic factor), ACTA2 (actin, alpha2, smooth skeletal muscle, aorta), ENSBTAG00000003151 (DNAJ heatshock protein family (Hsp40) member B11, GMPPB (GDP-mannose pyrophosphorylase B), GNPNAT1 (glucosamine-phosphate N-acetyltransferase 1). Three genes were down-regulated in the low RFI bulls. These included; LIMS2 (LIM zinc nger domain containing 2), GSTT1 (glutathione S-transferase theta 1), DBP (D-box binding PAR bzip transcription factor). HSPA5, CRELD2, HYOU1, ACTA2, DBP, MANF, and GSTT1 have previously been implicated in variation in feed e ciency [22-24, 34, 36, 45, 47-50] and are potential candidate biomarkers for this complex trait. Furthermore, HSPA5, CRELD2, HYOU1 and GSTT1 were all signi cantly associated and MANF tended towards a signi cant association with RFI phenotype further implicating the importance of these key genes to RFI, not only to the bulls used in the current study but to other cohorts of cattle divergent for RFI status. Gene ontology enrichment using DAVID revealed one enriched cluster in the up-regulated DEG from the low RFI bulls, which included signalling functions, disulphide bonding and peptide signalling. Ingenuity pathway analysis identi ed Aldosterone Signalling in Epithelial cells, GDP-mannose Biosynthesis, NRF2 mediated Oxidative Stress Response and EIF2 Signalling as the top canonical pathways related to RFI in bull hepatic tissue. The results of the GO enrichment and IPA analysis indicate that oxidative response, protein mediation and cell signalling in the liver are likely to be processes that are in uencing variation in feed e ciency.
The aldosterone signalling pathway was identi ed by IPA as the top canonical pathway due to the upregulation of the genes HSPA5 and DNAJB11 in the hepatic tissue of low RFI bulls. Aldosterone is secreted by the adrenal glands and has a major role in electrolyte and uid homeostasis. The aldosterone signalling pathway has been implicated in feed e ciency previously in the spleen of ine cient animals [31] and, interestingly, a GWAS study with cattle identi ed this pathway to be associated with variation in feed conversion ratio [51]. The protein encoded by HSPA5 is a member of the HSP70 family and as this protein interacts with many endoplasmic reticulum (ER) proteins it is likely to be important in monitoring protein transport through the cell [52]. DNAJB11 as a member of the DNA-J family of proteins is involved in the correct folding of proteins [53]. Speci cally DNAJB11 is involved in protein processing and metabolism of proteins and serves as a co-chaperone for HSPA5 in the ER [50]. One of the many functions of this family of proteins is to stabilize new proteins by ensuring correct folding or by helping refold proteins that have already been damaged by cell stress. Both DNAJB11 and HSPA5 have been implicated with feed e ciency previously [24,31,45,54]. In agreement with the current study, other work from our group [55] observed that animals undergoing compensatory growth with concomitant improvements in feed e ciency have a higher hepatic transcript abundance of DNAJB11 and HSPA5. Similarly, an up-regulation of HSPA5 was observed in the hepatic tissue of low RFI animals by Paradis et al. [34]. Taken together these results are indicative of a greater capacity in controlling cellular function and organisation as well as protein metabolism in more feed e cient animals.
NRF2 mediated oxidative stress response was also observed to be an overrepresented pathway in the hepatic tissue of high RFI bulls in the current study due to the up-regulation of the genes GSTT1 and DBP. NRF2 is a member of the cap 'n' collar basic region leucine zipper (cnc bZip) group of transcription factors [56]. This transcription factor is ubiquitously expressed in tissues but is only activated in response to a range of oxidative and electrophilic stimuli including ROS, antioxidants, glucose induced oxidative damage, heavy metals, and certain disease processes [56][57][58]. This canonical pathway has been associated with feed e ciency in a number of studies [29,49,50,54,59]. GSTT1 is a member of the glutathione S-transferase family and is involved in metabolism of xenobiotics and in catalysing reactions between the antioxidant glutathione and a host of potentially toxic compounds, highlighting it as an important homeostatic molecule [60]. The glutathione S-transferase family has previously been implicated in feed e ciency in various species [29,49,50,54,59]. In agreement with the current study Chen et al. [50] and Lindholm-Perry et al. [29] observed an up-regulation of these genes in feed ine cient cattle. DBP is a Protein Coding gene and amongst its cited functions is activation of circadian gene expression. Gene ontology annotations related to this gene include transcription factor activity, sequencespeci c DNA binding and transcriptional activator activity and RNA polymerase II core promoter proximal region sequence-speci c binding [61]. Additionally DBP was also implicated as contributing to variation in RFI phenotype through the network analysis reported by Weber et al. [23]. Similarly it has been shown that mice with increased FE have a lower expression of DBP [62] which is in agreement with the present study. Taken together, our work and that of the aforementioned authors, suggest that less e cient animals are exhibiting an increased oxidative stress, re ected in their increased anti-oxidation activities.
The differential expression and signi cant association of HSPA5, CRELD2, HYOU1, GSTT1 and MANF with RFI in the present study are also noteworthy as all six genes have been previously observed as differentially expressed in relation to feed e ciency in cattle [22,23,45,50,63]. While the biological signi cance of these genes in relation to feed e ciency remains unclear, due to their consistent presence in the literature, they should not be ruled out as potential biomarkers for this trait.

Conclusion
The present study contributes to the published knowledge base regarding the transcriptomic regulation of variation in feed e ciency. Our work, in combination with that of others as previously mentioned, highlights common genes underpinning feed e ciency in cattle, as measured by RFI, regardless of breed or genetic background. RNA-seq analysis is an exploratory approach that provides new hypotheses to be further investigated by other complementary approaches including global proteomics and ultimately, potential variation in the genes identi ed in this study may provide a basis for the selection of candidate biomarkers for the RFI trait and, following appropriate validation, contribute to genomic selection breeding programmes to improve feed e ciency in beef cattle. However, the key message from this work highlights the inconsistency in gene expression pro les across genders of the same genetic merit. While this inconsistency may be explained by the differences in physiological maturity of the two genders, it indicates that extensive further investigation is required before biomarker selection for RFI can be adopted.

Animal Model
The animals used in this study were derived from a purebred herd of Simmental cattle located at Teagasc Grange (Dunsany, Co. Meath, Irealnd) originally established to examine various aspects of the biological control of the RFI trait and which has been well characterised to date in the published literature [13][14][15][16][17][18][19]. In order to generate animals for the current study, the 20 highest (ine cient) and 20 lowest (e cient) ranking cows (n=40 in total) on phenotypic RFI were subjected to a multiple ovulation and embryo recovery programme and breeding using arti cial insemination to bulls with high and low estimated breeding values for RFI (see Crowley et al. [20]), respectively. There was no crossover of sires used across both cohorts of donor females. Following embryo recovery, all cows used were returned to the commercial herd within Teagasc Grange (Dunsany, Co. Meath, Ireland). Resultant embryos were transferred to crossbred beef heifer recipients (n=40), which were purchased from commercial herds speci cally for use in this study. In order to overcome potential confounding effects, animals were all purebred from the same genetic pool and were reared under the same conditions as a contemporary group from birth. All animals used were healthy, received regular veterinary interventions as required and were supervised on a daily basis by trained animal care staff. Pregnant heifers were managed under standard protocols and following calving were allowed to suckle their calves for a period of up to 7 days.
Following this, all heifers used were returned to the commercial herd within Teagasc Grange (Dunsany, Co. Meath, Ireland). In order to standardise rearing, calves (n=40) were then abruptly weaned and were subsequently reared on an electronic calf feeder. Brie y calves were offered milk replacer (MR) (Blossom Easymix; Volac, Co. Cavan, Ireland) and concentrate in pelleted form using an electronic feeding system (Vario; Foster-Tecknik, Engen, Germany), which recorded all feed-related events including intake of both

Pathway analysis
Biological processes, cellular components and molecular functions that were associated with the DEGs were identi ed using the GO enrichment tool of DAVID (Database for Annotation, Visualization and Integrated Discovery) Bioinformatics Resources v6.8 [65]. Fisher's exact test was used to determine the enrichment of the GO terms. Annotation clusters and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analyses were assessed using the o cial gene symbols of the DEGs and Bos Taurus as the species where possible. P-values associated with each annotation term inside each cluster are Fisher Exact/EASE Scores. To examine the molecular functions and genetic networks; the RNAseq data was further interrogated using IPA (Ingenuity Systems, Redwood City, CA; http://www.ingenuity.com), a web based software application that enables identi cation of over-represented biological mechanisms, pathways and functions most relevant to experimental datasets or genes of interest [55,66]. Data were imported in a exible format using the gene symbol as the identi er.