Effect of early calf-hood nutrition on the transcriptomic profile of subcutaneous adipose tissue in Holstein-Friesian bulls

Animal model

Tissue collection

The calves were euthanized at a mean age (±S.D.) of 126 (±1.1) days, using an intravenous overdose of sodium pentobarbitone. Death was confirmed by lack of ocular response. Subcutaneous adipose samples were obtained from the flank of the carcass following slaughter. The epidermis and dermis layers were cut to reveal the subcutaneous adipose. The tissue samples were washed with Dulbecco’s PBS and any blood vessels were removed. The tissue was subsampled and the subsamples were either (i) snap-frozen in liquid nitrogen and subsequently stored at − 80 °C for long term storage pending further analysis or (ii) fixed in 10% neutral buffered formalin (Sigma-Aldrich, Wicklow, Ireland) and embedded in paraffin in accordance to the standard procedures. Sections (5 μm thick) were dehydrated in ascending concentrations grades of alcohol, followed by clearing with xylene and were stained for gross anatomy using haematoxylin and eosin (H&E) stain.

Adipocyte cell number and diameter

Adipocyte cell number and cell size were analysed using Aperio ImageScope (v12; Leica Biosystems, Wetzlar, Germany) from the H&E stain slides. Five square sections (0.25 mm²) were electronically generated and randomly assigned on each slide and every adipocyte within these square sections was counted manually. The diameter of 200 adipocytes per animal, selected at random, within these square sections was measured at the widest point for each cell, using the ruler tool from the ImageScope software. The evaluator was blind to the treatments.

Statistical analysis on adipocyte cell size and number

Data for adipocyte cell size and number and bodyweight were analysed using the procedures of Statistical Analysis Software (SAS version 9.3, Cary, NC, USA) and were tested for normality (UNIVARIATE procedure). Residuals for this data were found to be normally distributed and were statistically analysed using ANOVA (MIXED procedure). Animal was the experimental unit and was included as a random effect. Sampling time (week of weight recording) was included in the statistical models as a repeated measure for weights. All results are presented as mean ± s.e.m.

RNA isolation and purification

Total RNA was extracted using RNeasy Lipid Tissue Mini kit from 100 mg of subcutaneous adipose (Qiagen, Manchester, UK). RNA was purified using the RNA Clean & Concentrator kit (Zymo Research, Irvine, CA, USA). The quantity of the RNA isolated was determined by measuring the absorbance at 260 nm using a Nanodrop spectrophotometer ND-1000 (Nanodrop Technologies, Wilmington, DE, USA). RNA quality was assessed on the Agilent Bioanalyzer 2100 (Agilent Technologies Ireland Ltd., Dublin, Ireland) using the RNA 6000 Nano Lab Chip kit (Agilent Technologies Ireland Ltd., Dublin, Ireland). Samples had an RNA Integrity Number (RIN) of mean (±SD) of 7.6 (±0.77).

RNAseq library preparation and sequencing

Twenty cDNA libraries were prepared from high quality RNA using an Illumina TruSeq RNA Sample Preparation kit v2 following the manufacturer’s instructions (Illumina, San Diego, CA, USA). For each sample, 1 μg of total RNA was used for cDNA preparation. All libraries were validated on the Agilent Bioanalyzer 2100 using the DNA 1000 Nano Lab Chip kit (Agilent Technologies Ireland Ltd., Dublin, Ireland). Individual RNAseq libraries were pooled (10 libraries per lane) based on their respective sample-specific-6 bp adaptors and sequenced at 100 bp/sequence single-end reads using an Illumina HiSeq 2500 sequencer.

RNAseq data analyses

Raw sequence reads were checked for quality using FASTQC software (version 0.10.0) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Input reads were then aligned to the bovine reference genome (UMD3.1) using Spliced Transcripts Alignment to a Reference (STAR) aligner. Reads were only aligned if they contained fewer than two mismatches with the reference genome and uniquely mapped to the reference genome. The software package HTSeq (v0.5.4p5) was used to calculate the number of sequence reads overlapping all protein coding genes from the ENSEMBLv74 annotation of the bovine genome. The number of read counts mapping to each annotated gene from HTSeq was then collated into a single file and used for calculation of subsequent differential gene expression. Only uniquely mapped reads were used for subsequent differential gene expression analysis.. Reads were filtered out if they contained fewer than five counts per million in nine samples out of the total sample number from subsequent analysis. The R (v 3.01) Bioconductor package EdgeR (v3.2.4), was applied to identify statistically significant differentially expressed genes (DEGs).

Pathway analysis

Biological pathways were identified using GOSeq software (v1.14.0) [17] and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway annotations [18] were used to identify biological pathways that were over-expressed among DEGs. GOSeq is an application for executing gene ontology analysis on RNAseq data while accounting for biases [17]. The online tool BioMart (www.ensembl.org/biomart/martview) was used to convert the bovine gene IDs, extracted from GOSeq, into human orthologs. These human orthologs were inputted into KEGG (http://www.genome.jp/kegg/pathway.html). To examine the molecular functions and biological pathways, the RNAseq data was also analysed using Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Redwood City, CA; https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/).

Animal performance and adipose tissue histology

Average daily bodyweight gain from the start to the end of the study was 1.08 ± 0.03 kg and 0.57 ± 0.03 kg for high and low plane of nutrition groups, respectively and this translated into a bodyweight difference of 53.8 kg at slaughter between the bulls of the high compared to the low treatment group (160.9 ± 3.98 kg versus 107.1 ± 3.19 kg, respectively; P < 0.001). Adipocyte cell number (291 ± 47.2 cells/1.25 mm² versus 7.5 ± 4.0 cells/1.25 mm², P < 0.001; Fig. 1) and cell diameter (36.9 ± 2.89 μm versus 3.7 ± 1.94 μm, P < 0.001) were both greater in the high compared with the low plane of nutrition animals, respectively.

Differential gene expression

The mean ± SD number of raw reads per sample was; 18,670,029 ± 4,093,759. Approximately 89% of the reads were uniquely mapped to the reference bovine genome after alignment. A multi-dimensional scaling (MDS) plot was created in Edge-R which estimated the degree of similarity between samples from the high compared to the low nutritional treatments. The MDS plot showed a separation between the two treatments (Fig. 2). There were 745 differentially expressed genes (DEG) between the two nutritional treatments (P < 0.05; False Discovery Rate < 0.05; fold change > 2.0). One hundred and sixty four genes had increased expression, while 581 had decreased expression in the low plane of nutrition compared to the high plane of nutrition. The RNAseq unprocessed data have been deposited in NCBI’s Gene Expression Omnibus and area accessible via GEO series accession number GSE97674.

Pathway analysis

Lipogenesis and adipogenesis

Animals on the low plane of nutrition had a lower of expression of MLXIPL (MLX interacting protein-like, also known as ChREBP; − 1.375 log fold change), a transcription factor that regulates lipogenesis and adipogenesis, in comparison to the high plane of nutrition. Target genes linked with adipocyte differentiation were down regulated in the low treatment in comparison to the high treatment, such as DGAT2 (diacylglycerol acyltransferase 2; − 5.215 log fold change), FASN (fatty acid synthetase; − 5-184 log fold change), ACACA (acetyl-CoA carboxylase-α; − 3.432 log fold change), LPL (lipoprotein lipase; − 3.241 log fold change), PLIN1 (Perilipin 1; − 2.279 log fold change) and FABP4 (fatty acid-binding protein 4; − 2.273 log fold change), INSIG-1 (insulin-induced gene 1; − 3.216 log fold change) and GDF-10 (growth differentiation factor 10; − 2.014). Genes that act as inhibitors for adipogenesis such as GATA2 (GATA binding protein 2; 1.151 log fold change) and WNT2B (Wnt family member 2B; 1.672 log fold change) were up regulated in the low plane of nutrition compared to the high.

Lipogenesis and adipogenesis

It has been postulated that rapid growth during the post weaning period in humans can cause preadipocytes to enter terminal differentiation [19]. Some transcription factors have been highlighted as regulators of lipogenesis and adipogenesis at a transcriptional level, including sterol regulatory element binding factor 1 (SREBF1) and MLXIPL [20]. The calves on the low plane of nutrition in our study had lower expression of MLXIPL in comparison to their contemporaries on the high plane of nutrition, which was supported by the complementary histology data. It has been reported that MLXIPL expression levels increase during preadipocyte differentiation in humans, mice and rats; however, there was a dramatic decrease in MLXIPL expression when cell differentiation took place in the absence of insulin [21, 22]. It is well accepted that insulin is a metabolically active hormone responsible for anabolic processes promoting growth, development, and nutrient homeostasis [23, 24]. Indeed we recorded greater plasma concentrations of insulin in the high compared with the low treatment groups at the time of slaughter (data not presented), although no differences were observed at the mRNA level. Insulin-induced gene-1 (INSIG-1) has been reported to mediate sterol regulatory element binding protein-1c (SREBP-c) which is a transcription factor involved in the regulation of the effect of insulin on lipid metabolism and adipogensis [25]. Studies have indicated that INSIG-1 increases in the fat tissue of mice developing diet induced obesity and also in differentiating 3 T3-L1 preadipocytes [26]. This highlights the importance of insulin in lipogenesis.

In a recent study [27] the authors reported a transcriptomic analysis of adipogenesis in human adipocyte tissue and identified a variety of genes as general adipocyte markers such as FABP4, PLIN1, LPL and ADIPOQ. These genes had increased expression during adipogenic differentiation consistent with the findings of others [28]. Upregulation of FABP4, FASN, ACACA and DGAT2 mRNA, which are target genes of PPARγ, are required for fatty acid biosynthesis and adipocyte differentiation [19, 29]. FASN is critical in the process of synthesising long chain saturated fatty acids and its function is inhibited by fasting [30]. Down regulation of ACACA has been reported in calorie restriction in adipose and liver tissue in pigs [31] and muscle in cattle [32]. It has been hypothesised that leptin controls adipocyte size via influencing DGAT expression [33]. Eventhough, no difference in PPARγ expression was determined in our study, downstream genes affected by PPARγ were downregulated in tissue from the animals on the low plane of nutrition. There was no evidence of BMP (bone morphogenic protein), expression being affected by plane of nutrition. However, GDF10 which is linked to BMP3 was found to downregulated on the low plane of nutrition. As well as a down regulation in genes that stimulate adipogenesis, in our study we found an upregulation in a key gene found to inhibit adipogenesis in the low plane of nutrition compared to the high plane of nutrition. It has been reported that GATA2 is expressed in preadipocytes and is down-regulated during terminal differentiation [34]. Studies in mouse preadipocytes show that GATA2 inhibits adipogenesis and traps cells at preadipocyte stage by direct suppression of PPARγ [34]. In addition, the histological data showing much lower adiposity (number and size of adipocytes) and greater number of preadipocytes in the animals on the low plane of nutrition corroborates the transcriptomic data aforementioned.

Endocrine regulatory networks

Insulin plays a role in regulating the synthesis and secretion of leptin, as adipocytes cultured in the presence of insulin increased the synthesis and secretion of leptin [35, 36]. It is therefore unsurprising that IPA predicted, when applied to our data, the stimulation of insulin as well as the known upregulation of LGALS12 (Galectin 12), BCAT2 (Branched Chain Amino Acid Transaminase 2) and ADIPOQ (adiponection), which are all located upstream of LEP and likely resulted in LEP to be upregulated in the high treatment. Plasma leptin is positively correlated with body fat in growing cattle [37] where it has been shown to signal energy abundance by informing the various bodily systems, including the HPT axis, of the prevailing metabolic status and energy reserves via circulating leptin protein concentrations [12]. For example, in adult animals, studies which examined the adipose tissue transcriptome of pre-partum dairy cows have shown that cows with greater BCS had greater expression of LEP [38]. Consistent with this we found that calves on the high plane of nutrition in our study had > 4.5 log fold higher LEP expression compared with animals the low animals. Despite this, when plasma concentrations were assayed, we failed to find any difference in leptin protein concentrations between the two groups (English et al., unpublished). This is in agreement with studies in both young dairy and beef bulls [1, 2, 39]. The lack of differences in plasma leptin concentrations may be due to the minute quantity of adipose tissue laid down in calves at this stage of development and/or the relative insensitivity of protein immunoassays compared with transcript based molecular approaches.

Leptin does not directly affect GnRH neurons of the hypothalamus and knocking out of leptin receptors on the GnRH neurons does not apparently delay the onset of puberty nor indeed affect subsequent fertility [13]. However, leptin has been found to indirectly influence the timing of puberty via the regulation of the hypothalamic Kiss1 neurons’ stimulation of GnRH [10]. Leptin can also bypass the hypothalamus and act on the pituitary gland and the testes; with leptin modestly stimulating gonadotropin secretion at a pituitary level and both direct stimulatory and inhibitory actions of leptin been reported in the gonads [40]. At high leptin concentrations, it has been reported to inhibit testosterone secretion in vitro in adult rats but to have no effect in pubertal rats [41]. This is possibly due to the fact that pubertal testes produce more 5α-reduced androgens than testosterone [42]. In research carried out by Byrne et al. (2016) bulls which had experienced dietary restriction prior to 6 months of age and an increase in plane of nutrition post 6 months of age leading to an increase in leptin concentration, were older at puberty than their contemporaries offered a high plane of nutrition; highlighting the importance of early life nutrition on age at puberty.

Adiponectin functions in the regulation of lipid and glucose metabolism, insulin sensitivity as well as inflammation [43–45]. Its concentration in blood has been shown to be negatively correlated to adipocyte size [46]. Food intake inhibits ADIPOQ mRNA expression in WAT and therefore, inhibits adiponectin concentration in serum [47]. Globular adiponectin protein has been found to inhibit the secretion of GnRH in GT1-7 cells derived from mouse hypothalamic GnRH neurons via the mediation of adenosine monophosphate-activated protein kinase (AMPK; Wen et al., 2008). However, high doses of adiponectin in MA-10 mouse Leydig cells has been reported to advance progesterone production through an increase in StAR (Steroidogenic Acute Regulatory Protein) and the CYP11A1 (Cytochrome P450 Family 11 Subfamily A Member 1) steroidogenesis enzyme, suggesting that adiponectin could stimulate testosterone production from the Leydig cells [48]. It has also been shown that adiponectin reduces LH secretion directly from the gonadotropes in the anterior pituitary via AMPK [49]. Pre-pubertal children with increased body mass index (BMI), though not classified as obese, had lower adiponectin concentration in subcutaneous adipose [50]. Similarly, postpartum cows, in negative energy balance had decreased ADIPOR1 and ADIPOR2 abundance in subcutaneous adipose compared to antepartum cows [51] which is in agreement with our findings with calves where we found ADIPOQ and its receptor ADIPOR2 were downregulated in the low compared to the high treatment groups. It has been reported that ADIPOR2 null mice demonstrated seminiferous tubules with aspermia but had normal testosterone concentration [52]. It has also been found that the abundance of ADIPOQ, ADIPOR1 and ADIPOR2 is greater in high fertility compared to that of medium and low fertility Holstein bulls [53]. There was no difference in serum concentrations of adiponectin between the two groups in the current study (English et al., unpublished), which like leptin may be due to the overall low adipose deposition in young bull calves. Using a very similar animal model and the same assays as those employed here Byrne et al. (2016); similarly observed in leptin and adiponectin concentrations during the same development as reported here. However, from 8 months on age onwards, bull calves offered a high plane of nutrition had greater leptin compared to their contemporaries on a moderate plane of nutrition (ADG c.0.65 kg/day).

Branched chain amino acid metabolism

Branched chain amino acids (BCAAs), leucine, isoleucine and valine are essential amino acids and, unlike other amino acids which are catabolized by the liver, BCAAs are catabolized mainly in the muscle, adipose, kidney and the brain; as the liver does not contain the branched amino acid transferase enzyme (BCAT2), necessary for BCAA catabolism [54]. The isoleucine/valine degradation pathway enters the TCA cycle either directly or via an acetyl derived intermediate and in this study we found that genes involved in the valine degradation pathway encoding BCAT2, 2-oxoisovalerate dehydrogenase (BCKDHB and DLD), 2-methylacyl-CoA dehydrogenase (ACADSB), enoyl-CoA hydratase (ECHS1, EHHADH, HADHA and HADHB), 3-hydroxyisobutyryl-CoA hydrolase (HIBCH), methylmalonate-semialdehyde dehydrogenase (ALDH6A1) and (S)-3-amino-2-methylpropionate transaminase (ABAT) were all down regulated in the low in contrast to the high treatment. Similarly, within the isoleucine degradation pathway, genes encoding for branched-chain-amino-acid transaminase (BCAT2), branched-chain α-keto acid dehydrogenase complex (DLD), 2-methylacyl-CoA dehydrogenase (ACADSB), enoyl-CoA hydratase (ECHS1, EHHADH, HADHA and HADHB), 3-hydroxy-2-methybutyryl-CoA dehydrogenase (HSD17B10) and acetyl-CoA C-acetyltransferase (ACAT1) were down regulated in the low compared to the high treatment. Expression of genes involved in branched chain amino acids metabolism were upregulated during 3 T3-L1 adipocyte differentiation [55, 56] suggesting that BCAA are involved in fatty acid synthesis [57]. It has also been shown that BCAA, especially leucine, systemic concentrations increase leptin secretion in vitro in murine [58] and decreases food intake in vitro in rats [59]. Systemic concentrations of branched chain amino acid are elevated in response to over nutrition in children and adolescents [60]. Therefore, it is reasonable to suggest that the evidence for a decline in BCAA degradation may be a result of the moderate restriction in dietary allowance experienced by our low plane of nutrition calves. The BCAA catabolic by-products mentioned in this study are intended for integration into the TCA cycle and therefore, any down regulation of the BCAA degradation pathway would cause a decrease in mitochondrial respiration in existing adipocytes [61].

Energy production and mitochondrial dysfunction