An examination of skeletal muscle and hepatic tissue transcriptomes from beef cattle divergent for residual feed intake

The selection of cattle with enhanced feed efficiency is of importance with regard to reducing feed costs in the beef industry. Global transcriptome profiling was undertaken on liver and skeletal muscle biopsies from Simmental heifers and bulls divergent for residual feed intake (RFI), a widely acknowledged feed efficiency phenotype, in order to identify genes that may be associated with this trait. We identified 5 genes (adj. p < 0.1) to be differentially expressed in skeletal muscle between high and low RFI heifers with all transcripts involved in oxidative phosphorylation and mitochondrial homeostasis. A total of 11 genes (adj. p < 0. 1) were differentially expressed in liver tissue between high and low RFI bulls with differentially expressed genes related to amino and nucleotide metabolism as well as endoplasmic reticulum protein processing. No genes were identified as differentially expressed in either heifer liver or bull muscle analyses. Results from this study show that the molecular control of RFI in young cattle is modified according to gender, which may be attributable to differences in physiological maturity between heifers and bulls of the same age. Despite this we have highlighted a number of genes that may hold potential as molecular biomarkers for RFI cattle.

Global agriculture is currently faced with the ambitious challenge of feeding a rapidly increasing global population, expected to peak at 9.2 billion by 2050 1 . This necessary increase in agricultural outputs must also be achieved within the current confines of arable land availability, thus it is essential that animal production systems become more efficient 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 importance 2 . Feed efficiency (dietary nutrient utilisation) in beef cattle is a trait of major economic importance 3 . Indeed, there is significant 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 efficiency it is possible to reduce feed intake in cattle while still maintaining growth and skeletal muscle gain and ultimately contribute to beef production profitability and sustainability.
Residual feed intake (RFI), first described by Koch et al. 7 is one such measure of feed efficiency, which can be defined 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, feed-efficient animals which have a low-RFI value and consume less than expected as well as feed inefficient animals which have a high-RFI and consume more than expected value may be identified 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 4 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-efficient, low-RFI animals and to proliferate their genetics through animal 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 goal for the identification of molecular biomarkers 20 . This not only provides an attractive alternative to direct measurement of dietary intake on large numbers of animals 6 , but it also allows for www.nature.com/scientificreports/ a better understanding of the biological mechanisms underlying RFI which is important 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 profiling, 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 20 . 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 21 . Furthermore, of these 5 genes only one, SCD, was differentially expressed in only one other molecular based RFI evaluation 22 . 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, particularly 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 23 . Thus, our efforts were focused towards examining the transcriptional alterations of these tissues between beef heifers and bulls divergent for high and low-RFI.

Results
Animal performance. Details of live weight gain, feed intake, and animal performance are presented in RNAseq read alignment and differential gene expression. RNA sequencing data was successfully generated for all samples with approximately 22.4 million sequences per sample available across both tissues.
On average approximately 85% of reads aligned to the bovine reference genome (UMD3.1) across both tissues analysed. The total number of genes expressed in each comparison was as follows: bull liver: 12,203 genes; bull muscle: 11,469 genes; heifer muscle: 10,766 genes and; heifer liver: 12,107 genes. A total of 5 genes (adj. p < 0.1) were identified as differentially expressed in skeletal muscle tissue between high and low RFI heifers with all 5 transcripts being up-regulated in the low RFI phenotype (Table 2). 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 were highly correlated with each other. A total of 11 genes (adj. p < 0.1) were identified as differentially expressed in hepatic tissue between high and low RFI bulls (Table 3) 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. Correlation analyses between DEGs and RFI value for each heifer and bull group, identified mostly significant associations. For example, DEGs from the heifer skeletal muscle analysis were all significantly negatively associated (p < 0.05) with RFI, with the exception of cytb, which tended towards significance (p = 0.0622). Similarly within the liver analysis of bulls divergent for RFI, 7 of the 11 DEGs were significantly associated with RFI (p < 0.05), with DBP, MANF and GMPPB tending towards significance (p < 0.1). ACTA2 expression was not significantly 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 4. None of the genes detected in heifer hepatic tissue and skeletal muscle tissue from bulls were found to be differentially expressed.
Differentially expressed genes and pathway analysis. Functional analysis of DEGs in hepatic tissue of bulls revealed amino sugar and nucleotide sugar metabolism as well as processes associated with the endoplasmic reticulum to be affected by divergence in RFI (Table 5). Functional analysis of heifer skeletal muscle DEGs indicated a role for oxidative phosphorylation and mitochondrial function towards divergence in RFI ( Table 6). All DEG were successfully mapped to a molecular or biological pathway and/or category in the Ingenuity Pathway analysis (IPA) database. DEGs 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 7. These included enriched pathways related to mitochondrial function and oxidative phosphorylation in heifers and Aldosterone Signalling, GDP-mannose Biosynthesis, Nuclear factor erythroid 2-related factor 2 (NRF2) mediated oxidative stress response, and Eukaryotic Initiation Factor 2 (EIF2) signalling in bulls.   www.nature.com/scientificreports/  www.nature.com/scientificreports/

Discussion
Residual feed intake is a prime feed-efficiency 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 20 . Furthermore, data from our own group provides evidence that RFI is a repeatable trait in beef cattle 11 . However, although well suited toward genetic propagation through genomic selection, insight into the underlying genes governing the trait are conflicting across studies 20,24 . This is undoubtedly due to the multifaceted nature of the trait as well as the influence 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 25 . Furthermore, RFI calculated within a specific population, as is the case in RFI transcriptional profiling studies, only reflects natural variation within that specific population, thus the level of divergence may be quite variable across different studies. In addition to these limitations, genomic selection models do not differentiate 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 specific 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, well-characterised 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 efficiency of an animal, with both organs being highly abundant in mitochondria 26 . 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 27 . 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 23 . We hypothesised that due to the metabolic importance of these organs, that variation in feed efficiency and energy expenditure (measured here using RFI) is likely to be reflected 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 significant 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 [28][29][30] . Therefore, the main observation of the current study is the inconsistent effect of RFI across (1) gender and (2) tissue for DEG profiles 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 there was some divergence in the actual physiological stage of development, as might be expected, between bulls and heifers at the time of sample collection. Similarly, in a targeted gene expression study, complementary to this current study, we observed a significant effect of gender on the expression of lipogenesis genes within the subcutaneous adipose tissue, this is despite there being no significant difference in subcutaneous fat measurements 31 . Typically, heifers display earlier physiological maturity when compared to bulls of the same age, thus the identification 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, it is possible that greater differences may have been apparent at the molecular level, had sample-size and sequencing depth been greater. Moreover, the low number of genes identified in the current study 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 influencing 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 www.nature.com/scientificreports/ evidence for a relationship between RFI and the aldosterone signalling pathway and the NRF2 mediated oxidative stress pathway in hepatic tissue of bulls. We identified 5 genes that were statistically 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 findings of other previously published RFI based studies [32][33][34] . 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 up-regulated in low RFI heifers compared to their high RFI counterparts. Additionally, all of these genes were significantly associated with RFI, with the exception of cytb, which only tended towards a significant association, further indicating 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 efficiency. 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 DEGs related to the complexes of ETC and RFI 12,35 . For example, increased levels of COXII (complex IV) and NADH dehydrogenase subunits have been shown to be associated with efficient animals at both the protein and transcriptome level 12,35-38 . However, although genes of the ETC have previously been reported in relation to variation in RFI, this is the first 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 specific 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. 21 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 type 21 .
Using gene ontology (GO) enrichment analysis, we identified important processes underlying feed efficiency variation in skeletal muscle of heifers. These included functions related to mitochondrial metabolism, in particular oxidative phosphorylation. Similarly, the top canonical pathways identified 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 efficiency across varying species 12,35,[39][40][41][42] . 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 43 . It seems likely that variation in mitochondrial function could contribute to variation in energy utilisation. Furthermore, it has been demonstrated that feed efficient animals exhibit greater capacity to modulate conditions of oxidative stress 44 . Feed efficient animals have been shown to have a higher activity of all enzymes of the ETC across multiple species including broilers and lambs 39,40,42 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 efficient heifers. These results complement previous research suggesting an association between decreased respiration capacity and increased ROS production in less efficient animals 40,45 . Moreover, the results of the present study are reinforced by the observation of Kong et al. 38 in which the mitoproteome was skewed towards high feed efficiency birds despite no difference in mitochondrial DNA between phenotypes, suggesting an increase in mitochondrial activity in the high feed efficiency phenotype 38 .
Within the hepatic tissue of bulls we identified 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 (glucosaminephosphate N-acetyltransferase 1). Three genes were down-regulated in the low RFI bulls. These included; LIMS2 (LIM zinc finger 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 efficiency [21][22][23]33,35,44,[46][47][48][49] and are potential candidate biomarkers for this complex trait. Furthermore, HSPA5, CRELD2, HYOU1 and GSTT1 were all significantly associated and MANF tended towards a significant 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 analysis of hepatic tissue in bulls identified terms including those related to amino sugar and nucleotide sugar metabolism as well as others related to endoplasmic reticulum protein processing as significantly enriched. Additionally, Ingenuity pathway analysis identified 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 processing and cell signalling in the liver are likely to be processes that are influencing variation in feed efficiency.
The aldosterone signalling pathway was identified by IPA as the top canonical pathway due to the up-regulation 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 fluid homeostasis. The aldosterone signalling pathway has been implicated in feed efficiency previously in the spleen of inefficient animals 30 and, interestingly, a GWAS analysis with cattle identified this pathway to be associated with variation in feed conversion ratio 50  www.nature.com/scientificreports/ 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 51 . DNAJB11 as a member of the DNA-J family of proteins is involved in the correct folding of proteins 52 . Specifically DNAJB11 is involved in protein processing and metabolism of proteins and serves as a co-chaperone for HSPA5 in the ER 49 .
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 cellular stress. Both DNAJB11 and HSPA5 have been implicated with feed efficiency previously 23,30,44,53 . In agreement with the current study, other work from our group 54 observed that animals undergoing compensatory growth with concomitant improvements in feed efficiency 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. 33 . Taken together these results are indicative of a greater capacity in controlling cellular function and organisation as well as protein metabolism in more feed efficient 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 55 . 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 [55][56][57] . This canonical pathway has previously been associated with feed efficiency in a number of studies 28,48,49,53,58 . 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 59 . The glutathione S-transferase family has previously been implicated in feed efficiency in various species 28,48,49,53,58 . In agreement with the current study, Chen et al. 49 and Lindholm-Perry et al. 28 observed an up-regulation of these genes in feed inefficient 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, sequence-specific DNA binding and transcriptional activator activity and RNA polymerase II core promoter proximal region sequence-specific binding 60 . Additionally, DBP was also implicated as contributing to variation in RFI phenotype through the network analysis reported by Weber et al. 22 . Similarly, it has been shown that mice with increased FE have a lower expression of DBP 61 which is in agreement with the present study. Taken together, our work and that of the aforementioned authors, suggest that less efficient animals are exhibiting an increased oxidative stress, reflected in their increased anti-oxidation activities.
The differential expression and significant 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 efficiency in cattle 21,22,44,49,62 . While the biological significance of these genes in relation to feed efficiency 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 efficiency. Our work, in combination with that of others as previously mentioned, highlights common genes underpinning feed efficiency 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 identified 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 efficiency in beef cattle. However, the key message from this work highlights the inconsistency in gene expression profiles 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.

Methods
All procedures involving animals in this study were conducted under an experimental licence from the Health Products Regulatory Authority in accordance with the cruelty to Animals Act 1876 and the European Communities (Amendment of Cruelty to Animals Act 1876) Regulation 2002 and 2005.
Animal model. The animals used in this study were derived from a purebred herd of Simmental cattle 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] . Briefly, purebred Simmental beef cattle derived from a herd phenotypically ranked on RFI were used for the current study 31 . Based on RFI phenotype the highest (n = 20) and lowest (n = 20) ranking cows and heifers were bred to pedigree Simmental sires with estimated breeding values for high and low RFI (see Crowley et al. 4 ), respectively, through artificial insemination and multiple ovulation/embryo recovery technologies. There was no crossover of sires used across both cohorts of donor females. Resultant embryos were transferred to crossbred beef heifer recipients and the resulting calves (bulls, n = 16; heifers n = 18) were used for the current study. Pregnant heifers were managed under standard protocols and following calving were allowed to suckle their calves for a period of up to 7 days. In order to standardise rearing, calves were then abruptly weaned and were subsequently reared on an electronic calf feeder. Briefly 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 RNAseq data analysis. FASTQC software (v0.11.5) was used to check the quality of the raw sequencing reads. Input reads were then aligned to the bovine reference genome (UMD3.1) using STAR (v2.5.1). HTSeq (v0.6.1p2) software was used to calculate the number of sequenced fragments overlapping all protein-coding genes from the ENSEMBLv88 annotation of the bovine genome. The number of counts of reads mapping to each annotated gene from HTSeq was then collated into a single file and used for subsequent differential gene expression. The R (v3.4) Bioconductor package, EdgeR (v3.5) 64 , which uses a negative binomial distribution model to account for both biological and technical variability, was applied to identify statistically significant DEGs. The analysis was undertaken using moderated tagwise dispersion. An adjusted p value < 0.1 (Benjamini-Hochberg adjustment) was applied as a threshold to call genes with differential expression levels. The R (v3.4) Corrplot package 65 was used to visualise the correlation matrices pertaining to the DEGs passing multiple correction for each comparison. Gene expression results for DEGs identified through transcriptional profiling were also www.nature.com/scientificreports/ correlated with RFI value, to determine potential associations between the DEGs and RFI. Correlations were undertaken using the CORR procedure in SAS 9.1 (SAS Inst. INC., Cary, NC).
Gene ontology and pathway analyses. Biological processes, cellular components and molecular functions that were associated with the DEGs were identified using the g:Profiler GO enrichment tool 66 . The gene name and associated p-value were used as input for gene ontology and pathway analyses. Fisher's exact test was used to determine the enrichment of the GO terms. 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. ingen uity. com), a web based software application that enables identification of over-represented biological mechanisms, pathways and functions most relevant to experimental datasets or genes of interest 67 . Data were imported in a flexible format using the gene symbol as the identifier.
Ethics declarations. All procedures involving animals in this study were conducted under an experimental licence from the Health Products Regulatory Authority in accordance with the cruelty to Animals Act 1876 and the European Communities (Amendment of Cruelty to Animals Act 1876) Regulation 2002 and 2005. All procedures described in this study were carried out in compliance with the ARRIVE guidelines.

Data availability
The datasets generated and analysed during the current study are available in the National Centre for Biotech- nology  www.nature.com/scientificreports/