RNA-seq analysis of bovine adipose tissue in heifers fed diets differing in energy and protein content

Adipose tissue is no longer considered a mere energy reserve, but a metabolically and hormonally active organ strongly associated with the regulation of whole-body metabolism. Knowledge of adipose metabolic regulatory function is of great importance in cattle management, as it affects the efficiency and manner with which an animal converts feedstuff to milk, meat and fat. However, the molecular mechanisms regulating metabolism in bovine adipose tissue are still not fully elucidated. The emergence of next-generation sequencing technologies has facilitated the analysis of metabolic function and regulation at the global gene expression level. The aim of this study was to investigate the effect of diets differing in protein and energy density level on gene expression in adipose tissue of growing replacement dairy heifers using next-generation RNA sequencing (RNAseq). Norwegian Red heifers were fed either a high- or low-protein concentrate (HP/LP) and a high- or low-energy roughage (HE/LE) diet from 3 months of age until confirmed pregnancy to give four treatments (viz, HPHE, HPLE, LPHE, LPLE) with different growth profiles. Subcutaneous adipose tissue sampled at 12 months of age was analyzed for gene expression differences using RNAseq. The largest difference in gene expression was found between LPHE and LPLE heifers, for which 1092 genes were significantly differentially expressed, representing an up-regulation of mitochondrial function, lipid, carbohydrate and amino acid metabolism as well as changes in the antioxidant system in adipose tissue of LPHE heifers. Differences between HPHE and HPLE heifers were much smaller, and dominated by genes representing NAD biosynthesis, as was the significantly differentially expressed genes (DEG) common to both HE-LE contrasts. Differences between HP and LP groups within each energy treatment were minimal. This study emphasizes the importance of transcriptional regulation of adipose tissue energy metabolism, and identifies candidate genes for further studies on early-stage obesity and glucose load in dairy cattle.


Introduction
The global human population is steadily increasing, and with it so are the demands for efficient animal production systems [1,2]. Within the dairy industry, efficient meat and milk production is dependent on animals with a well-functioning metabolism and an efficient and purposeful nutrient allocation [3][4][5]. For growing cattle this implies a rapid growth rate without excess fattening, as the number of days spent to achieve productive size is proportional to total rearing costs [6], while excessive fattening has a negative impact on both feed efficiency, subsequent health and carcass value [5,7]. High growth rates can be achieved by both breeding and feeding strategies. Growth is also affected by management and environmental conditions, although to a lesser degree [8,9]. As the capacity for rapid genetic progress at a farm level is limited and management and environment yield smaller effects, nutritional management is the main strategy applied in practice. However, nutritionally induced high growth rates are associated with increased fatness [10]. There is a plethora of heifer studies reporting the effects of feeding regimes during the rearing period on growth rate and characteristics, as well as on subsequent milk yield in dairy heifers [11,12]. However, the molecular mechanisms by which these effects are exerted are still not fully elucidated.
An important role of adipose tissue in this regard has emerged, as adipose is no longer considered merely a passive energy storage tissue, but is a transcriptionally active endocrine organ with an important role in whole-body metabolism and energy balance [13,14]. When compared to monogastrics this is especially true for ruminants, in which adipose tissue is the main site for both lipid synthesis, storage and mobilization [15]. Increased knowledge of how transcriptomic regulation in adipose tissue is affected by nutrition in dairy heifers may contribute to the development of both feeding and breeding strategies to produce more efficient animals. Several studies on adipose tissue from cattle fed differently have been performed in the recent years, elucidating the connection between plane of nutrition and the expression of genes associated with lipogenesis, adipogenesis and lipolysis [16][17][18]. However, most gene expression studies on adipose tissue in cattle to date have been performed using microarrays and targeted gene analysis [14][15][16]. Recent increased availability and reduced costs of next generation sequencing techniques have led to more gene expression studies being performed by RNAsequencing (RNAseq) technology. Compared to microarrays, RNA-seq is more sensitive and has a much wider dynamic range. It analyzes the whole transcriptome in an unbiased fashion, and is therefore not limited by predetermined target genes, array setups or a need for existing genome annotations [19].
Data in the published literature provides evidence of an effect of nutrition on adipose expression of genes, but until now such information have only been available through microarrays and targeted gene expression analyses performed on adult cows [16][17][18]. Hence, the objective of this study was to investigate the effects of diets differing in protein and energy density on the full adipose transcriptome of growing dairy heifers, and to identify gene networks related to critical metabolic pathways which can be targeted in further studies to enhance replacement heifer rearing strategies and production efficiency in the Norwegian Red cow.

Ethics statement
All experimental procedures involving animals were approved by the Norwegian Animal Research Authority (FOTS ID 2955, Reference number: 2010/203231). Biopsies were harvested under epidural anesthesia, and all animals were given 3 mg/kg BW of Ketoprofen intramuscularly (Comforion vet., Orion Pharma Animal Health) to prevent pain or inflammation at the biopsy site.

Animal model
This study was part of a large experiment designed to examine the effects of diets with differing energy and protein contents on growth characteristics, fat deposition and subsequent milk yield in Norwegian Red heifers, and from this to determine the optimal growth curve in the pre-and post-pubertal periods for a replacement dairy heifer. Eighty heifers from the dairy herd at the Norwegian University of Life Sciences (years 2010 and 2011) were assigned either to a high (HE) or low (LE) energy group, fed according to a BW gain of 800-1000 or 600-750 g/day from three months of age to confirmed pregnancy, respectively. Each energy group was further subdivided into two protein groups, i.e. low (LP) or high (HP) to give four dietary treatments with 20 animals in each group, viz low-protein high-energy (LPHE), high-protein high-energy (HPHE), low-protein low-energy (LPLE) and high-protein low-energy (HPLE). The university herd consists of two genetic groups, one group with high milk-yield (HMY) and one group with low occurrence of clinical mastitis (LCM) [20]. As this genetic difference could affect several variables of interest to our experiment, it was taken into account when allocating heifers to dietary treatments to generate groups which were balanced according to genetic group. Apart from this consideration, heifers were randomly allocated to treatment group. Body weights [21] were measured every second week, and body condition scores (BCS) and height at withers (WH) recorded every month for all experimental animals. The Norwegian 1-5 scale BCS system is based on Edmonson et al. [22] but calibrated for Norwegian Red [23].
Heifer feeding and management. All heifer calves were fed the same diet from birth to three mo of age. From three months of age until confirmed pregnancy, heifers were housed in a tie-stall barn and fed one of the four experimental diets. Fresh feed was offered twice daily and individual feed intake was recorded four days a week. HE groups were fed grass silage ad libitum whereas LE groups were fed the same silage mixed with 10 to 40% wheat straw (on a dry matter basis) given as restricted rations. All heifers were daily fed 1 kg of either of two custom-made concentrates, differing in protein content. Energy density of the diets was adjusted with the roughage quality. Diet protein content was adjusted with both roughage quality and type of concentrate. Table 1 gives the average nutrient content of roughage and custom made concentrates. Table 2 gives average daily intake of dry matter, energy, and protein for the different treatment groups.

Adipose tissue biopsy sampling
Adipose tissue biopsies were collected at 12 months of age. The heifers were given a low epidural anesthesia of 3-5 ml of a local anesthetic (Lidokel-Adrenalin vet, Kela Laboratoria NV) prior to sampling. Subcutaneous adipose tissue biopsies were harvested lateral to the tail base. Each biopsy consisted of 0.5-1.5 g of adipose tissue. Biopsies were immediately snap frozen in liquid nitrogen, and subsequently stored at -80˚C.

RNA extraction and sequencing
Each adipose tissue biopsy was crushed on dry ice, and 100 mg was used for RNA extraction. RNA extraction was carried out using the Qiagen RNeasy lipid tissue mini kit (Qiagen). All samples were homogenized in Trizol using a PRO 200 handheld homogenizer (PRO Scientific Inc., Oxford, CT USA). Except for the homogenization step, the RNA extraction methodology was performed according to the manufacturer's protocol. The quantity and quality of RNA was determined using a Nanodrop spectrophotometer ND-1000 (Nanodrop Technologies, Wilmington, DE, USA) and the RNA 6000 Nano Lab Chip kit (Agilent Technologies Ireland Ltd., Dublin, Ireland), respectively. The best-quality RNA samples to give a balanced sample set with respect to genetic group and dietary treatment were selected for further RNA sequencing (n = 6 for LPHE, HPHE, LPLE and HPLE, respectively). RIN values for these samples varied between 7.8 and 9.2. To avoid DNA contamination, all RNA samples were subjected to a DNAse treatment after extraction (Turbo DNA-free, Ambion, Life Technologies) and a subsequent clean-up (RNA Clean and Concentrator -25, Zymo Research Corp.).
Library preparation was performed using the Illumina TruSeq RNA sample prep kit v2, according to the manufacturer's instructions (Illumina, San Diego, CA, USA). Sequencing was carried out on an Illumina HiSeq 2000 workstation, with 100 bp paired end reads, and the 48 samples distributed randomly over 12 lanes (Clinical Genomics, Toronto).

Statistics and bioinformatics
The fastq files containing the raw sequence reads were examined using the software FASTQC (v 0.11.2.). Adapters were trimmed from all sample fastq files using Cutadapt (v 1.3). Reads were simultaneously quality trimmed by removing reads shorter than 20 bp, or with Phred scores below 25. Following adapter and quality trimming, all files were checked again in FASTQC. Paired-end read files were aligned to the UMD3.1 bovine genome assembly [24], using Tophat 2.0.12 [25]. The number of reads per gene in each sample was counted using HTSeq-count (v 0.6.1) [26]. Statistical analysis of read counts were carried out in EdgeR (v 3.1.2), a Bioconductor software package run in the statistical software environment R (v 3.1.2) [27]. Differential gene expression was tested in EdgeR using tagwise dispersions. The statistical tests were corrected for multiple testing using the Benjamini-Hochberg method as implemented in EdgeR. A generalized linear model which included dietary treatment and genetic group was fitted to the data using the glmFit function [28]. The comparisons were performed Table 1. Average nutrient content of high-and low-energy (HE and LE) roughages and high and low protein (HP and LP) concentrates, all weighted by number of feed days during the experimental feeding period from 3 mo of age to confirmed pregnancy (g/kg of dry matter if not stated otherwise). as pairwise comparisons and the differences between contrasting treatments, i.e. LPHE-LPLE, HPHE-HPLE, LPHE-HPHE and LPLE-HPLE, were tested. Significantly differentially expressed genes (DEG) were called at a FDR (false discovery rate adjusted p-value as given in EdgeR) of 0.1, and retained for further analysis. Subsequently, the resulting lists of DEG from the low and high-energy HP/LP contrasts and the low and high-protein HE/LE contrasts were compared to identify any DEG common to both comparisons within an energy or protein level. This suggests this gene was affected by energy or protein level regardless of the level of the other dietary variable. The DEG were subsequently analyzed using Ingenuity 1 Pathway Analysis (IPA 1 , QIAGEN Redwood City, www.qiagen.com/ingenuity) to identify affected biological pathways and cellular functions represented by the DEG. An affected pathway is defined as a pathway containing more DEG than expected by chance, given the number of genes belonging to the pathway. Significance level for the pathway analyses was set at p < 0.05. Pathways or functions specifically pertaining to irrelevant diseases, species or tissues were omitted.

Animal performance
At 12 months, BW of heifers were significantly different (p < 0.05) between all dietary treatments, but less so between different protein treatments within an energy treatment than across energy treatments. Body condition scores were different between different energy treatments (p < 0.05), but not between protein treatments within an energy treatment. Height at withers were significantly different (p < 0.05) between all treatment groups, but the largest difference was found between different energy levels. Least-squares means of BW, BCS and WH at 12 months of age are shown in Table 3.

Differential gene expression
Sequencing yielded on average 57.2 million raw reads per sample. Alignment analysis resulted in an average overall read mapping rate of 95.1% and a concordant pair alignment rate of 89.6%. Genes with a minimum CPM (counts per million) of two in at least three samples were Table 2 considered to be expressed and were included in the differential expression analysis. Out of 26,740 gene transcripts in the UMD3.1 bovine genome assembly, 12,476 gene transcripts fulfilled this criterion. A multidimensional scaling plot visualizing the level of similarity between dietary treatments is presented in Fig 1. The greatest difference in adipose gene expression profiles was found between LPHE and LPLE heifers. This contrast yielded 1092 DEG, 40 times more than any other comparison. The number of DEG identified by contrasting treatments, energy or protein levels are shown in Table 4. Contrasts between the HP and LP groups within each energy level yielded only two DEG for the HE contrast, and one DEG for the LE contrast. Details on these DEG are presented in Table 5. The full lists of DEG between the LPHE-LPLE and the HPHE-HPLE contrasts are shown in S1 and S2 Tables, respectively. The RNAseq raw data are deposited and accessible through NCBI's Gene Expression Omnibus with the GEO Series accession number GSE79347 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE79347).

Pathway analysis
Differentially expressed genes were grouped in IPA according to their biological functions. The functional groups most highly represented by DEG found for LPHE-LPLE contrasts were nucleic acid metabolism, amino acid metabolism and small molecule biochemistry. Small molecule biochemistry was also the functional group most highly represented by the HPHE-HPLE DEG. Figs 2 and 3 display further cellular function differences found for DEG from the two LPHE-LPLE and HPHE-HPLE contrasts, respectively. Pathway analysis of the DEG between LPHE and LPLE heifers resulted in 99 significantly affected pathways (p < 0.05). The main pathways affected by the difference in LPHE and LPLE feeding regimes pertained to mitochondrial functions (mitochondrial dysfunction, oxidative phosphorylation and TCA cycle), amino acids (valine, isoleucine, leucine, tryptophan, arginine, proline and cysteine metabolism and degradation) and carbohydrate metabolism (pyruvate, propanoate, butanoate and ketone body metabolism, glycolysis and gluconeogenesis), antioxidant defense systems and fatty acid biosynthesis and degradation. The most highly affected pathways for this contrast are shown in Fig 4. Pathway analysis of the 24 DEG from the HPHE-HPLE comparison yielded the IPA pathways shown in Fig 5. Significantly affected pathways for this contrast were mainly related to glycosaminoglycan and NAD biosynthesis. NAD biosynthesis was also the main pathway represented by the 17 DEG common to both HE-LE contrasts. DEG from the two HP-LP contrasts were not subjected to pathway analysis due to low numbers.

Gene expression differences between energy treatments
The massive difference in DEG numbers identified for the examined contrasts suggests a profound interaction between energy and protein level in the diet on adipose tissue metabolism, and that this interaction effect is in fact larger than the effects of protein or energy independently. The differences in gene expression between LPHE and LPLE heifers illustrates the importance of transcriptional regulation in metabolism, as several of the main affected pathways such as oxidative phosphorylation, TCA cycle, amino acid and carbohydrate metabolism clearly reflect changes in energy metabolism which are expected to differ in animals fed highor low-energy diets. The DEG between HPHE and HPLE heifers suggest the main differences  RNA-seq of adipose tissue in heifers between these two groups was related to connective tissue structure (upregulated glycosaminoglycan synthesis by xylosyltransferase 1, XYLT1), and NAD synthesis by quinolinate phosphoribosyltransferase (QPRT). The low number of DEG in this contrast also affects the number of DEG common to both HE-LE contrasts, and leaves NAD synthesis as the pathway mainly affected by energy level regardless of protein level in the diet. This is an interesting finding, although not surprising, as nicotinamide adenine dinucleotide (NAD) is a coenzyme well known to take part in a plethora of cellular oxidation-reduction reactions, which constitute the basis of cellular energetics [29]. Based on the results of other gene expression studies, we expected certain genes to be differentially expressed in both energy contrasts. These genes included FASN [18,30], SCD [18,30], ELOVL6 [30], ACACA [18,30] and OB [31]. We would like to emphasize that even though these "classical" and expected DEG from the LPHE-LPLE contrast, did not obtain an FDR < 0.1 in the HPHE-HPLE contrast, 991 of the 1092 DEG from the LP comparison displayed the same directional change in the HP contrast. Thus, differences were mainly the same, only smaller, and given a bigger sample set, the results for the two contrasts would probably be more similar. For comparative purposes, we also performed  an overall comparison of all HE animals against all LE animals within the same fitted model as described above. With an FDR cutoff value of 0.05 (because of the increased power of this test), the resulting DEG list was strikingly similar to the DEG list for the LPHE-LPLE contrast. In the following 4 sub-sections the main differences in adipose gene expression between LPHE and LPLE heifers are discussed. Mitochondrial function. The tricarboxylic acid (TCA) cycle and the electron transport chain performing oxidative phosphorylation take place in the mitochondria. The TCA cycle is central to many anabolic and catabolic pathways in a mammalian cell as its main substrate, acetyl-CoA is produced from the initial breakdown of carbohydrates, protein and lipids. It also provides substrates for the synthesis of fatty acids, sterols, amino acids and glucose. The cycle consists of eight steps. Genes associated with all eight steps were up-regulated in LPHE heifers which indicates increased mitochondrial citrate production in LPHE heifers. Citrate is the main fatty acid precursor originating from oxidative metabolism. Citrate destined for fatty acid synthesis is transported out of the mitochondria by a citrate transport protein encoded by SLC25A1, which was up-regulated in LPHE heifers. In the cytosol, citrate is converted to acetyl-CoA by ATP citrate lyase (ACLY). In ruminants, acetate (from ruminal fermentation) has been shown to be the main carbon source for fatty acid synthesis [32]. Citrate lyase has traditionally been regarded as virtually inactive in ruminant adipose tissue, because of the low levels RNA-seq of adipose tissue in heifers of glucose available for metabolism in peripheral tissues [32]. However, ACLY was up-regulated in LPHE compared to LPLE heifers. Recent findings, including this study show that the ACLY gene displays expression differences in the adipose tissue of cattle with different energy states [18,30]. This suggests that citrate lyase does indeed play a part in ruminant metabolism and metabolic regulation, and may indicate an increased incorporation of glucose into fatty acids in animals with a positive energy balance. Genes associated with all of the 5 complexes in the electron transport chain were up-regulated in the LPHE heifers, indicating an increase in electron transport chain activity and energy (ATP) production. This makes sense as the role of adipose tissue as a metabolic excess energy endpoint requires it to deal with all forms of energy-yielding nutrients, transforming them in to lipids or to the energy required to produce these lipids. The genes associated with changes in mitochondrial function between LPHE and LPLE heifers, are presented in Fig 6. Fatty acid biosynthesis and degradation. Differentially expressed genes associated with fatty acid synthesis and degradation are shown in Table 6. The up-regulation of genes associated with fatty acid synthesis (ACACA, SCD, FASN and ELOVL6), fatty acid activation (ACSL4, SLC27A3, SLC27A6), and triacylglycerol (TAG) synthesis (AGPAT2, DGAT2, ELOVL6, LPCAT4, PLPP3) in the LPHE heifers is a reflection of their increased energy intake, leading to a higher production of TAG from citrate and other precursors than in the LPLE heifers. The activity of delta-9-desaturase is important in ruminant adipose tissue compared to monogastrics. This is because of the extensive microbial biohydrogenation of dietary unsaturated fatty acids occurring in the rumen, leaving almost exclusively saturated fatty acids (SFA) to be transported from viscera to the peripheral tissues in ruminants [33]. However, the high melting point of SFA makes them unsuitable as the sole source of fatty acids incorporated into cell membranes and adipose tissue. The high delta-9-desaturase activity of ruminant adipose tissue alleviates this effect, and makes oleic acid the most abundant FA in ruminant tissues [34].
Amino acid metabolism. Differentially expressed genes pertaining to amino acid metabolism are shown in Table 7. Genes associated with the metabolism and degradation of several amino acids (isoleucine, leucine, phenylalanine, valine, cysteine, lysine, alanine, tryptophan and aspartate) were up-regulated in the LPHE heifers. This is likely a general response, which would be present also in e.g. muscle because HE animals grow much faster than the LE heifers and thus have both an increasingly higher daily crude protein intake as well as an increased metabolism and degradation of amino acids. Carbohydrate metabolism. Differentially expressed genes pertaining to carbohydrate metabolism are shown in Table 8. Carbohydrates are energy-yielding nutrients, and several genes associated with carbohydrate degradation such as glycolytic reactions and mitochondrial transport and breakdown of pyruvate were up-regulated in the LPHE heifers. The fact that glycolysis-associated genes were DE between the two energy levels, suggests that the LPHE heifers had sufficient glucose available to direct some of it towards adipose tissue energy storage. This is not considered a common fate for glucose in ruminants which under most physiological conditions are dependent on de novo gluconeogenesis in the liver, as most of the glucose in the diet will be degraded by rumen microbes [32]. However, more recent studies have shown that hepatic de novo gluconeogenesis is highly prioritized in cows, and proceeds at a relatively high rate even in the absence of lactational demands or in the presence of available dietary glucose [35]. Thus, glucose directed towards adipose energy storage in HE heifers may originate from hepatic synthesis. The pentose-phosphate pathway is mainly associated with anabolic metabolism. It's main functions are to provide ribose-5-phosphate for nucleic acid biosynthesis and reducing equivalents (NADPH) for the biosynthesis of fatty acids and steroids, as well as for the reduction of reactive oxygen species (ROS) [36]. The up-regulation of genes coding for pentose-phosphate pathway enzymes in HE heifers is in concordance with the up-regulation of genes associated with fatty acid synthesis and glutathione redox reactions (Tables 6 and 9, respectively). The up-regulation of genes associated with propionate, butyrate and ketone body metabolism and transport in LPHE heifers indicate an increased flux/transport of propionate and 3-hydroxybutyrate from blood into adipocytes. This is possible because the LPHE heifers had an increasingly higher energy intake per day than the LPLE heifers, and as a result microbial fermentation in ruminants causes the energy in feed to enter the animal's bloodstream mainly as volatile fatty acids (VFA), not as glucose or fatty acids [37]. The difference in energy intake between HPHE and HPLE heifers was comparable to the LP energy contrast, however there was no significant up-regulation of genes pertaining to VFA transport between these two groups. A possible explanation is that the increased fraction of protein-rich feedstuffs in the experimental concentrate given to the HP groups, was substituted by barley, which is rich in starch, in the low-protein concentrate fed to the LP groups. This led to a difference in daily starch intake of approximately 130 grams between protein treatments. Starch is rapidly degraded into VFA (mainly propionate) in the rumen. The increased amount in the LP diets combined with the energy surplus in the HE roughage may have resulted in a sufficiently increased VFA flux from the rumen to the blood to evoke an up-regulation of 11 genes in LPHE heifers, enabling their adipose tissue to deal with the energy excess in the form of VFA and direct this energy towards fat deposition. Antioxidant functions. Differentially expressed genes pertaining to cellular antioxidant functions are shown in Table 9. Virtually any metabolic reaction involving oxygen may yield reactive oxygen species (ROS) as byproducts. Consequently, both mitochondria, endoplasmic reticulum (ER) and peroxisomes are major sites of ROS production in the cell [38]. In controlled quantities, these substances have a function in cell signaling and regulation [38]. However, accumulation of ROS will lead to a state of oxidative stress, causing uncontrolled oxidation of DNA and proteins in the cell and cellular dysfunction [39]. Therefore, it is imperative that cells have well-functioning systems to regulate and control the level of ROS at all times. The glutathione system plays an important role in the cellular defense against oxidative damage and harmful substances. Glutathione is synthetized by glutamate cysteine ligase (GCL, encoded by GCLC) and glutathione synthetase (GSS, encoded by GSS) from glycine and cysteine. It functions as an electron donor for reactive oxygen species (ROS), itself being dimerized to glutathione disulfide (GSSG) in the reaction. Glutathione disulfide may be recycled back to 2 GSH molecules by glutathione reductase, using NADPH as an H+ donor. In concordance with studies performed on mouse models, there was an up-regulation of GCLC and GSS and a down-regulation of the glutathione peroxidase gene GPX7 in LPHE heifers [40,41]. When examined in isolation, a down-regulation of GPX7 indicates a decreased level of redox reactions producing GSSG and is associated with increases in oxidative stress [40], inflammatory status and insulin resistance [41]. However, as shown in Fig 6, a whole range of other antioxidant system genes were upregulated in LPHE heifers. Among these were several genes coding for glutathione-S-transferases (GSTT1, GSTT3, GSTM1, MGST1 and MGST3). These enzymes Table 6. Fatty acid synthesis and degradation-associated genes differentially expressed between low-protein, high-energy heifers and low-protein, low-energy heifers. represent another part of the glutathione system. They are multifunctional, and may act as peroxidases as well as catalyzing the conjugation of glutathione with unwanted substrates such as epoxides and lipid hydroperoxides, making them eligible for transport out of the cell [42]. The shift towards a conjugate-forming oxidative defense rather than one relying solely on redox reactions involving GSH/GSSG also indicates an increased cellular export of glutathione in HE heifers, as these conjugates are excreted and catabolized extracellularly [43]. An increased export of glutathione conjugates out of the cell may explain the need for increased production of GSH and the up-regulation of GCLC and GSS. An increased cellular excretion of harmful substances could also be expected from the upregulation of EPHX2. This gene codes for soluble epoxide hydrolase (SEH), a peroxisomal enzyme playing a role in cellular detoxification by converting lipid epoxides to trans-hydro-diols, which are eligible for conjugation and subsequent excretion [44]. The genes SOD2 and CAT were also up-regulated in HE heifers. Superoxide dismutase (SOD) converts superoxides to H 2 O 2 , which may subsequently be reduced by either GPX or catalase. If regarded in isolation, the up-regulation of these two genes probably would represent a deterioration of the antioxidant defense in the HE heifers, as catalase exclusively binds H 2 O 2 , and at a low affinity, while GPX has the ability to reduce several other organic peroxides as well. For this reason, catalase has been regarded as a second line antioxidant defense, taking effect when peroxide levels become supraphysiological [45].

Symbol
The combination of DEG representing mitochondrial energy metabolism pathways such as the TCA cycle and oxidative phosphorylation as well as the antioxidant system, is what makes the mitochondrial dysfunction pathway so highly ranked in the LPHE-LPLE contrast. During the last decade, increasing evidence has highlighted the associations between nutritional overload, mitochondrial changes in redox and energy state, insulin resistance and impaired adipocyte function [39][40][41]45]. Although our HE heifers were clinically healthy, it is possible that these heifers, at 12 months of age and a BCS of almost 4 were on the verge of a harmful degree of fatness, and thus displayed a gene expression profile with similarities to that of morbidly obese individuals of other species. If this is the case, we may have revealed candidate genes for the early detection of obesity. Among these are genes coding for glutathione-S-transferases (GSTT1, GSTT3, GSTM1, MGST1 and MGST3), SOD2 and CAT.

Expression differences between protein treatments
The differences between the protein treatments within each energy group were minimal. Only CRYM, coding for μ-crystallin, was DE for the HPLE-LPLE contrast, and was downregulated in the HPLE heifers. The expression of CRYM is negatively associated with serum glucose levels, and positively with insulin sensitivity [46]. Therefore, it is interesting that this gene was DE in heifers fed different protein levels, but not different energy levels. Depending on passage rate and heat treatment, the ruminal digestibility of ground barley (present in varying amounts in our experimental concentrates) varies from 80-90% [47,48], and the main product is propionate. Thus, only 10-20% of ingested starch may pass undegraded from the rumen to the intestine and can be absorbed into the blood as glucose. Given the higher starch intake of LPLE heifers, one would expect a higher blood glucose and a decreased CRYM expression in these animals, but it is possible that the association of both CRYM expression and blood glucose levels with insulin sensitivity may actually cause the opposite effect if insulin sensitivity is sufficiently high. If a difference in daily starch intake has indeed resulted in the observed up-regulation of CRYM, this indicates that certain parts of the metabolic regulation in heifers are sensitive to small changes in protein and carbohydrates available through feed, even on isocaloric diets, and that CRYM may be a sensitive indicator of glucose load in cattle. For the HPHE-LPHE contrast only DSP and ACTA1 were DE, for which any notable association to dietary protein or starch are as of yet unknown. Thus, we have not found any genes to be differentially expressed as an effect of different feed protein levels regardless of energy level. This emphasizes further the dominant effect of total ration composition above the single effects of dietary energy or protein on the adipose tissue gene expression of growing heifers.

Conclusions
Dietary energy content has a massive impact on the bovine adipose tissue transcriptome regulating several aspects of energy metabolism. The size of the effect is modulated by dietary protein content. The different protein levels applied in this study had only minor effects on adipose gene expression per se. Compared to heifers with a BCS of 3.63, heifers fed to a BCS of 3.95 at 12 mo of age displayed several transcriptional responses which have been associated with increased oxidative stress in other species. This study is a first step towards the identification of biomarkers for early-stage obesity (GSTT1, GSTT3, GSTM1, MGST1, MGST3, SOD2, CAT) and glucose load (CRYM) in cattle. However, certain biomarker identification will require further validation studies such as functional assays, lipidomics and proteomics. The current study underpins the importance of adipose transcriptional regulation of metabolism in growing cattle and highlights how total ration characteristics have a larger impact on adipose tissue transcriptomes than the single effects of protein and energy.
Supporting information S1