# Transcriptomic changes in liver of young bulls caused by diets low in mineral and protein contents and supplemented with n-3 fatty acids and conjugated linoleic acid

Sara Pegolo, Alessio Cecchinato, Núria Mach, Massimiliano Babbucci, Marianna Pauletto, Luca Bargelloni, Stefano Schiavon, Giovanni Bittante

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sara Pegolo, Alessio Cecchinato, Núria Mach, Massimiliano Babbucci, Marianna Pauletto, et al. Transcriptomic changes in liver of young bulls caused by diets low in mineral and protein contents and supplemented with n-3 fatty acids and conjugated linoleic acid. PLoS ONE, 2016, 11 (12), pp.1-24. 10.1371/journal.pone.0167747. [hal-01602348](https://hal.science/hal-01602348)

HAL Id: hal-01602348

https://hal.science/hal-01602348

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

CrossMark

click for updates

RESEARCH ARTICLE

# Transcriptomic Changes in Liver of Young Bulls Caused by Diets Low in Mineral and Protein Contents and Supplemented with n-3 Fatty Acids and Conjugated Linoleic Acid

Sara Pegolo¹*, Alessio Cecchinato¹, Núria Mach², Massimiliano Babbucci³, Marianna Pauletto³, Luca Bargelloni³, Stefano Schiavon¹, Giovanni Bittante¹

1 Department of Agronomy, Food, Natural Resources, Animals and Environment (DAFNAE), University of Padova, Legnaro, Padova, Italy, 2 Animal Genetics and Integrative Biology unit (GABI), INRA, AgroParisTech, Université Paris-Saclay, Jouy-en-Josas, France, 3 Department of Comparative Biomedicine and Food Science, University of Padova, Legnaro, Padova, Italy

* sara.pegolo@unipd.it

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Pegolo S, Cecchinato A, Mach N, Babbucci M, Pauletto M, Bargelloni L, et al. (2016) Transcriptomic Changes in Liver of Young Bulls Caused by Diets Low in Mineral and Protein Contents and Supplemented with n-3 Fatty Acids and Conjugated Linoleic Acid. PLoS ONE 11(12): e0167747. doi:10.1371/journal.pone.0167747

**Editor:** Robert S. Weiss, Cornell University, UNITED STATES

**Received:** July 21, 2016

**Accepted:** November 19, 2016

**Published:** December 8, 2016

**Copyright:** © 2016 Pegolo et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** Raw Illumina sequencing data have been deposited in GenBank (SRA) (PRJNA303144, BioProject SRP066556, study accession). All other relevant data are within the paper and its Supporting Information files.

**Funding:** This research was supported by the University of Padova (Ricerca Scientifica fondi quota EX 60% -649 60A08-7943/14).

**Competing Interests:** The authors have declared that no competing interests exist.

## Abstract

The aim of the present study was to identify transcriptional modifications and regulatory networks accounting for physiological and metabolic responses to specific nutrients in the liver of young Belgian Blue × Holstein bulls using RNA-sequencing. A larger trial has been carried out in which animals were fed with different diets: 1) a conventional diet; 2) a low-protein/ low-mineral diet (low-impact diet) and 3) a diet enriched in n-3 fatty acids (FAs), conjugated linoleic acid (CLA) and vitamin E (nutraceutical diet). The initial hypothesis was that the administration of low-impact and nutraceutical diets might influence the transcriptional profiles in bovine liver and the resultant nutrient fluxes, which are essential for optimal liver function and nutrient interconversion. Results showed that the nutraceutical diet significantly reduced subcutaneous fat covering *in vivo* and liver pH. Dietary treatments did not affect overall liver fat content, but significantly modified the liver profile of 33 FA traits (out of the total 89 identified by gas-chromatography). In bulls fed nutraceutical diet, the percentage of n-3 and CLA FAs increased around 2.5-fold compared with the other diets, whereas the ratio of n6/n3 decreased 2.5-fold. Liver transcriptomic analyses revealed a total of 198 differentially expressed genes (DEGs) when comparing low-impact, nutraceutical and conventional diets, with the nutraceutical diet showing the greatest effects on liver transcriptome. Functional analyses using ClueGo and Ingenuity Pathway Analysis evidenced that DEGs in bovine liver were variously involved in energy reserve metabolic process, glutathione metabolism, and carbohydrate and lipid metabolism. Modifications in feeding strategies affected key transcription factors regulating the expression of several genes involved in fatty acid metabolism, e.g. insulin-induced gene 1, insulin receptor substrate 2, and RAR-related orphan receptor C. This study provides noteworthy insights into the molecular changes occurring as a result of nutrient variation in diets (aimed at reducing the environmental impact and improving human health) and broadens our understanding of the relationship between nutrients variation and phenotypic effects.
Introduction

Growing attention has been recently paid to sustainability in the livestock sector, the function-
ality of animals, and product quality according to the EU Framework Program Horizon 2020
[1]. Sustainable beef production depends on a balanced combination of various factors, such
as the reduction of negative environmental impacts, animal health, nutritional efficiency, effi-
cient growth, and a desirable carcass composition with a high percentage of lean tissue and a
low fat content. The Nitrate Directive (91/676/EC), which regards the protection of water qual-
ity across Europe, together with the increasing cost of protein sources have brought pressure
on farmers to reduce the N content in livestock diets thereby reducing N excretion and emis-
sion [2–3]. Minerals, which are often given to animals in excess of their requirements [4], can
also contribute to environmental pollution (especially P, Cu and Zn).

Aside from environmental issues, a growing awareness of the association between diet and human health has increased the demand for foods of high nutritional quality. There is now recognition of the importance of fatty acids (FAs) of the n-3 series and conjugated linoleic acid (CLA), the predominant isomers being 18:2c9,t11 (9-CLA, the rumenic acid, RA) and 18:2t10, c12 (10-CLA), for their potential beneficial effects on the functioning of the immune, nervous, and cardiovascular systems in humans [5], as well as for reproductive performances and carcass traits in ruminants [6]. Dietary supplementation with n-3 FAs or CLA has been found to affect FA composition in ruminant tissues [7–9]. Specific effects have been ascribed to CLA mixtures, where two isomers generally dominate, i.e. 18:2c9,t11 and 18:2t10,c12. The former isomer has anti-carcinogenic properties [10], whereas the latter can increase the lean to fat ratio in growing animals [11–12]. Adaptation to low protein and minerals, and to high n-3 FAs and CLA ingestion is influenced by the transcriptional and translational regulation of genes that encode the proteins controlling these processes. To our knowledge, studies are available about the effects of these FAs on the expression of target genes in ruminants [13–15] but no information is available on the effects of protein and mineral deficiency.

Transcriptomic approaches have been applied for genome-wide screening in blood and tis-
sue biopsies aimed at better understanding the gene expression networks that control nutrient
utilization. In beef cattle, few studies adopted RNA-sequencing methodology aiming to
improve the understanding of the supplement-driven physiological and metabolic mecha-
nisms that regulate the global tissue response to nutrients [16–18].

In ruminants, *de novo* fatty acid synthesis predominantly takes place in the adipose tissue
but also the liver has a specific role in lipid metabolism [19]. In particular, ruminant liver is
involved in the uptake, oxidation and metabolism of non-esterified fatty acids (NEFAs), syn-
thesis of cholesterol and phospholipids, and synthesis and secretion of lipoproteins. Polyunsat-
urated fatty acids (PUFAs), after rumen biohydrogenation, could be further metabolized in the
liver and their subsequent deposition in muscle may be affected [20]. Thereby, modification of
muscle fatty acid composition aiming to enhance the nutritional value of meat should take
into consideration liver metabolism. We hypothesized that identification of the liver genes spe-
cifically regulated by n-3 and CLA, as well as the variations in dietary protein and mineral con-
tents could reveal unique biomarkers of physiological adaptation to diet supplementation and
could also provide significant insights into the molecular control of this response and the
resulting phenotypes. Based on the phenotypic results of a large experimental trial, involving
males and females of Belgian Blue × Holstein and Belgian Blue × Brown Swiss, we selected the
12 Belgian Blue × Holstein bulls (because of their greater importance for the beef industry in
Europe) to perform a comparative transcriptome analysis of liver in bulls fed with a conven-
tional diet (Conv), a low-protein/low-mineral diet (low-impact diet, LowI) or a low-protein/
low-mineral diet supplemented with n-3 FAs, CLA, and vitamin E (nutraceutical diet, Nutr)
using RNA-sequencing. The specific aims of this work were: (i) to describe the changes in liver gene expression profiles and phenotypic variations when comparing a Conv with a LowI or a Nutr diet, (ii) to identify the biological processes in liver regulated by variations in diet composition, (iii) to define the relationships between the expression profiles of regulated genes and phenotypic traits such as FA profiles.

# Results

## Phenotypic analyses

**Growth performances, carcass and liver traits.** No significant differences among treatments were found for body weight (BW), average daily gain (ADG), and *in vivo* muscularity score (Table 1).

In contrast, significantly lower values of *in vivo* fat covering were found with Nutr compared with Conv and LowI groups (*P* < 0.05). None of the carcass traits considered (carcass weight and yield, muscularity and fatness scores) differed significantly among treatments, even though muscularity and fatness scores were nominally lower in Nutr compared with

<table><caption>Table 1. Growth performances, carcass traits, liver characteristics of young bulls fed conventional, low-impact and nutraceutical diets.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Diet<sup>1</sup></th><th rowspan="2">SEM</th><th rowspan="2">P-values<sup>2</sup></th></tr><tr><th>Conv</th><th>Lowl</th><th>Nutr</th></tr></thead><tbody><tr><td><strong>Growth performances</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Final body weight, kg</td><td>684</td><td>671</td><td>662</td><td>17.23</td><td>0.893</td></tr><tr><td>Average daily gain, kg</td><td>1.54</td><td>1.42</td><td>1.38</td><td>0.03</td><td>0.366</td></tr><tr><td>Muscularity score<sup>3</sup></td><td>4.22</td><td>4.17</td><td>4.25</td><td>0.07</td><td>0.916</td></tr><tr><td>Fatness score<sup>3</sup></td><td>1.89<sup>a</sup></td><td>1.83<sup>a</sup></td><td>1.41<sup>b</sup></td><td>0.06</td><td>0.042</td></tr><tr><td><strong>Carcass traits</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Carcass weight, kg</td><td>396</td><td>388</td><td>382</td><td>9.49</td><td>0.580</td></tr><tr><td>Carcass yield, g/kg</td><td>0.580</td><td>0.578</td><td>0.576</td><td>0.002</td><td>0.619</td></tr><tr><td>Muscularity score<sup>4</sup></td><td>4.00</td><td>4.00</td><td>3.75</td><td>0.12</td><td>0.788</td></tr><tr><td>Fatness score<sup>4</sup></td><td>2.00</td><td>2.08</td><td>1.83</td><td>0.11</td><td>0.641</td></tr><tr><td><strong>Liver</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Weight, kg</td><td>6.53</td><td>6.36</td><td>6.90</td><td>0.15</td><td>0.750</td></tr><tr><td>Dry matter, %</td><td>29.23</td><td>29.61</td><td>29.45</td><td>0.15</td><td>0.515</td></tr><tr><td>Fat, %</td><td>3.06</td><td>2.45</td><td>2.38</td><td>0.22</td><td>0.358</td></tr><tr><td>Protein, %</td><td>19.85</td><td>19.23</td><td>19.67</td><td>0.12</td><td>0.149</td></tr><tr><td>Ash, %</td><td>1.53</td><td>2.09</td><td>1.77</td><td>0.07</td><td>0.066</td></tr><tr><td>pH</td><td>6.13<sup>a</sup></td><td>6.14<sup>a</sup></td><td>6.08<sup>b</sup></td><td>0.01</td><td>0.044</td></tr></tbody></table>

SEM: Mean standard error of the means.

¹Mean values are reported for each diet.

²P-values after Kruskal Wallis non-parametric test

³In accordance with Schiavon et al. [21], body conformation was linearly scored in vivo from S+ (all muscle profiles extremely convex; exceptional muscle development) to P- (all muscle profiles concave to very concave; poor muscle development) taking into account the profiles of shoulders, loins, rump, thighs, and buttocks (S+ = 6.33; P- = 0.66), and fat covering was linearly scored *in vivo* in 5 classes (1 = very lean; 5 = very fat) using a combined visual and palpation approach taking into account the presence and thickness of subcutaneous fat depots at the base of the tail, ribs, and shoulders.

⁴SEUROP scoring system for carcass muscularity from S+ (all muscle profiles extremely convex; exceptional muscle development) to P- (all muscle profiles concave to very concave; poor muscle development) taking into account the profiles of shoulders, loins, rump, thighs, and buttocks (S+ = 6.33; P- = 0.66), and for carcass fat covering (1 = very lean; 5 = very fat).

Letters indicate significant differences (*P* < 0.05; Bonferroni post-hoc test).

doi:10.1371/journal.pone.0167747.t001
Conv and LowI treatments. No differences in liver weight among treatments were observed. Average liver pH was significantly lower when animals were fed the Nutr diet than when fed the LowI or Conv diet.

**Liver fatty acid and mineral profile.** Nineteen FAs out of the total 89 identified by GC showed significant differences (*P* < 0.05) among dietary treatments (Table 2).

Only some minor individual saturated FAs were affected (*P* < 0.05) by dietary treatments (17:0*iso* and 24:0). The proportions of some monounsaturated FAs (MUFAs) were affected by the Nutr diet. In particular, oleic acid (18:1*c9*) was lower in Nutr young bulls compared with Conv (-29%) and LowI (-23%) treatments (*P* = 0.050). A notable reduction in 18:1*cl1*, and some medium-chain MUFAs (14:1 isomers) was also observed in Nutr vs Conv (*P* = 0.042 and *P* = 0.034, respectively), while increased proportions of 18:1*cl2* (*P* = 0.034) and other 18:1 isomers were detected (*P* = 0.018).

The relative content of some PUFAs was affected by the Nutr diet (Table 2). For instance, significant differences (*P* = 0.038) were detected for the linoleic acid (18:2 *c9*,*c12*) content, which represented almost half of all PUFAs in the liver of LowI bulls, and only 39% of Conv and 37% of Conv bulls. Compared with Conv, Nutr strongly increased the content of alpha-linolenic acid (18:3*c9*,*c12*,*c15*; +365%) and other minor 18:3 isomers (*P*<0.05). The CLA *c*,*t*,*t*,*c* isoforms and total CLA isomers increased about 2.5-fold in Nutr compared with Conv and LowI (*P* = 0.050). About individual CLA contents, large differences were found for 18:2*c9*,*c11*, 18:2*tl0*,*c12* and 18:2*cl1*,*t13* (*P*<0.05). The n-3 FAs proportion also increased 2.5-fold in Nutr, with a contemporaneous reduction in the n-6/n-3 ratio in Nutr compared with the other groups (*P* = 0.043). In addition, n-3 FAs eicosapentaenoic acid (EPA) and docopentaenoic acid (DPA) were much higher (*P* = 0.018) in Nutr compared with Conv and LowI. *Trans* 18:1 FAs significantly increased in the young bulls fed Nutr (1.8-fold, *P* = 0.043) compared with Conv.

The content of 27 minerals in liver was determined by inductive couple plasma-optical emission spectroscopy (ICP-OES), and of these 27 minerals, 11 (As, B, Ba, Co, Cr, Hg, Li, Ni, PB, Sb, and V) were below the limit of detection. Aluminum and Si were excluded due to a number of missing data. Liver concentrations were detected for 14 minerals: Ca, Cd, Cu, Fe, K, Mg, Mn, Mo, Na, P, S, Se, Sn, and Zn (Table 3), but no significant differences were observed among the 3 feeding groups.

Transcriptomic data analyses

**Sequencing and mapping.** RNA-sequencing data were obtained from 3 Conv bulls, 4 LowI bulls and 4 Nutr bulls. All samples passed quality control measures for raw and trimmed sequenced reads. An average of over 20 million reads were obtained per sample, with ~80% of the reads mapping to the bovine reference genome (see summary statistics in S1 Table). High correlations between read count tools for total and unique reads were obtained (R² > 0.99). A boxplot of gene counts after TMM normalization revealed no bias in count distribution across samples (S1 Fig). Before performing the differential gene expression analysis, further multivariate analyses were carried out on the normalized gene count data to identify any discrepancies within and between individuals. The MDS plot, a technique aiming to detect meaningful underlying dimensions that allow explaining observed dissimilarities (distances) between samples, was also applied to the 24,596 genes after the TMM normalization. Inspection of the MDS plot showed Conv and LowI samples to partially overlap while a clearer distinction emerged for Conv and LowI vs Nutr (Fig 1).

In addition, sample BC7_Conv appeared to deviate markedly from the other samples. The interclass PCA also revealed that bull BC7_Conv laid outside the 95% confidence interval of its
<table><caption>Table 2. Liver fatty acid (FA) contents¹ of young bulls fed conventional (Conv), low-impact (Lowl) and nutraceutical (Nutr) diets.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Diet<sup>2</sup></th><th rowspan="2">SEM</th><th rowspan="2">P-values<sup>4</sup></th></tr><tr><th>Conv</th><th>Lowl</th><th>Nutr<sup>3</sup></th><th>SEM</th></tr></thead><tbody><tr><td colspan="7"><strong>Individual FAs, g/100g FAs</strong></td></tr><tr><td>17:0<sup>iso</sup></td><td>0.336<sup>a</sup></td><td>0.291<sup>ab</sup></td><td>0.249<sup>b</sup></td><td>0.006</td><td>0.034</td></tr><tr><td>17:0<sup>anteiso</sup></td><td>0.696</td><td>0.481</td><td>0.400</td><td>0.016</td><td>0.051</td></tr><tr><td>24:0</td><td>0.215<sup>ab</sup></td><td>0.147<sup>b</sup></td><td>0.278<sup>a</sup></td><td>0.009</td><td>0.034</td></tr><tr><td>14:1 others<sup>5</sup></td><td>0.023<sup>ab</sup></td><td>0.034<sup>a</sup></td><td>0.004<sup>b</sup></td><td>0.002</td><td>0.042</td></tr><tr><td>15:1 sum<sup>5</sup></td><td>0.006</td><td>0.024</td><td>0.000</td><td>0.001</td><td>0.058</td></tr><tr><td>17:1 sum<sup>5</sup></td><td>0.261</td><td>0.219</td><td>0.158</td><td>0.006</td><td>0.051</td></tr><tr><td>18:1<sup>t12</sup></td><td>0.451</td><td>0.677</td><td>0.721</td><td>0.017</td><td>0.055</td></tr><tr><td>18:1<sup>c9</sup></td><td>8.887</td><td>8.208</td><td>6.348</td><td>0.209</td><td>0.050</td></tr><tr><td>18:1<sup>c11</sup></td><td>1.025<sup>a</sup></td><td>0.919<sup>ab</sup></td><td>0.682<sup>b</sup></td><td>0.011</td><td>0.034</td></tr><tr><td>18:1<sup>c12</sup></td><td>0.349<sup>b</sup></td><td>0.501<sup>ab</sup></td><td>0.585<sup>a</sup></td><td>0.012</td><td>0.043</td></tr><tr><td>18:1 others<sup>6</sup></td><td>0.400<sup>b</sup></td><td>0.682<sup>ab</sup></td><td>1.308<sup>a</sup></td><td>0.028</td><td>0.018</td></tr><tr><td>18:2<sup>c9</sup>,<sup>c12</sup></td><td>15.932<sup>b</sup></td><td>20.181<sup>a</sup></td><td>16.129<sup>b</sup></td><td>0.160</td><td>0.038</td></tr><tr><td>18:3<sup>c9</sup>,<sup>c12</sup>,<sup>c15</sup></td><td>0.620<sup>b</sup></td><td>0.759<sup>ab</sup></td><td>2.263<sup>a</sup></td><td>0.033</td><td>0.034</td></tr><tr><td>18:3 others<sup>5</sup></td><td>0.003<sup>b</sup></td><td>0.000<sup>b</sup></td><td>0.620<sup>a</sup></td><td>0.061</td><td>0.022</td></tr><tr><td>18:2<sup>c11</sup>,<sup>t13</sup> (CLA)</td><td>0.028<sup>b</sup></td><td>0.034<sup>b</sup></td><td>0.129<sup>a</sup></td><td>0.006</td><td>0.034</td></tr><tr><td>18:2<sup>c9</sup>,<sup>c11</sup> (CLA)</td><td>0.014<sup>b</sup></td><td>0.013<sup>b</sup></td><td>0.073<sup>a</sup></td><td>0.008</td><td>0.030</td></tr><tr><td>18:2<sup>t10</sup>,<sup>c12</sup> (CLA)</td><td>0.000<sup>b</sup></td><td>0.016<sup>b</sup></td><td>0.132<sup>a</sup></td><td>0.005</td><td>0.028</td></tr><tr><td>20:3<sup>c8</sup>,<sup>c11</sup>,<sup>c14</sup></td><td>4.784</td><td>2.740</td><td>3.104</td><td>0.147</td><td>0.064</td></tr><tr><td>20:4<sup>c8</sup>,<sup>c11</sup>,<sup>c14</sup>,<sup>c17</sup></td><td>0.315<sup>a</sup></td><td>0.212<sup>a</sup></td><td>0.882<sup>ab</sup></td><td>0.014</td><td>0.043</td></tr><tr><td>20:4 others<sup>5</sup></td><td>0.120</td><td>0.116</td><td>0.255</td><td>0.004</td><td>0.055</td></tr><tr><td>20:5<sup>c5</sup>,<sup>c8</sup>,<sup>c11</sup>,<sup>c14</sup>,<sup>c17</sup> (EPA)</td><td>0.332<sup>ab</sup></td><td>0.245<sup>b</sup></td><td>0.920<sup>a</sup></td><td>0.011</td><td>0.018</td></tr><tr><td>22:4<sup>c7</sup>,<sup>c10</sup>,<sup>c13</sup>,<sup>c16</sup></td><td>3.753<sup>a</sup></td><td>2.744<sup>a</sup></td><td>1.995<sup>ab</sup></td><td>0.083</td><td>0.018</td></tr><tr><td>22:5<sup>c4</sup>,<sup>c7</sup>,<sup>c10</sup>,<sup>c13</sup>,<sup>c16</sup></td><td>0.884</td><td>0.610</td><td>0.613</td><td>0.019</td><td>0.055</td></tr><tr><td>22:5<sup>c4</sup>,<sup>c7</sup>,<sup>c10</sup>,<sup>c13</sup>,<sup>c16</sup>,<sup>c19</sup> (DPA)</td><td>1.929<sup>ab</sup></td><td>1.445<sup>b</sup></td><td>3.734<sup>a</sup></td><td>0.033</td><td>0.018</td></tr><tr><td>20:2 to 20:5 others<sup>5</sup></td><td>0.065</td><td>0.017</td><td>0.082</td><td>0.006</td><td>0.065</td></tr><tr><td>22:6<sup>c4</sup>,<sup>c7</sup>,<sup>c10</sup>,<sup>c13</sup>,<sup>c16</sup>,<sup>c19</sup> (DHA)</td><td>0.322</td><td>0.292</td><td>0.617</td><td>0.016</td><td>0.050</td></tr><tr><td colspan="7"><strong>Groups of FA, g/100g FAs</strong></td></tr><tr><td>PUFA</td><td>40.767</td><td>40.832</td><td>44.142</td><td>0.256</td><td>0.057</td></tr><tr><td>CLA c,t/t,c</td><td>0.365</td><td>0.409</td><td>0.894</td><td>0.481</td><td>0.050</td></tr><tr><td>CLA total</td><td>0.445</td><td>0.474</td><td>1.053</td><td>0.024</td><td>0.050</td></tr><tr><td>Trans 18:1</td><td>1.623<sup>b</sup></td><td>2.418<sup>ab</sup></td><td>2.896<sup>a</sup></td><td>0.078</td><td>0.043</td></tr><tr><td>n-6</td><td>35.296</td><td>35.477</td><td>32.031</td><td>0.253</td><td>0.055</td></tr><tr><td>n-3</td><td>3.995<sup>ab</sup></td><td>3.923<sup>b</sup></td><td>9.696<sup>a</sup></td><td>0.116</td><td>0.043</td></tr><tr><td>n-6/n-3</td><td>8.855<sup>ab</sup></td><td>9.257<sup>a</sup></td><td>3.327<sup>b</sup></td><td>0.124</td><td>0.043</td></tr></tbody></table>

SEM: Mean standard error of the means; CLA: conjugated linoleic acid; PUFA: polyunsaturated fatty acids; EPA: eicosapentaenoic acid; DPA: docopentaenoic acid; DHA: docohexaenoic acid

Letters indicate significant differences (P < 0.05; Bonferroni post-hoc test)

¹Only the traits with significant differences (P < 0.05) or trends (P < 0.10) are reported

²Mean values are reported for each diet.

³One animal was not included in the FA analysis due to bad quality of the chromatography data.

⁴P-values after Kruskal Wallis non-parametric test.

⁵Sum of positional isomers

⁶18:1 others: sum of 18:1t4, 18:1t6, 18:1t8, 18:1t9, 18:1t10, 18:1t13 + t14, 18:1c13, 18:1t16, 18:1c14, 18:1c15, and 18:1c16.

doi:10.1371/journal.pone.0167747.t002
<table><caption>Table 3. Mineral contents of diets and liver of young bulls fed conventional (Conv, n = 3), low-impact (LowI, n = 4) and nutraceutical diets (Nutr, n = 4).</caption><thead><tr><th rowspan="2"></th><th colspan="3">Diet<fnref n="1" /></th><th rowspan="2">SEM<fnref n="2" /></th><th rowspan="2">P-values<fnref n="2" /></th></tr><tr><th>Conv</th><th>LowI</th><th>Nutr</th></tr></thead><tbody><tr><td colspan="6"><strong>Micro-element, mg/kg</strong></td></tr><tr><td>Cd</td><td>0.065</td><td>0.070</td><td>0.066</td><td>0.003</td><td>0.47</td></tr><tr><td>Cu</td><td>91.2</td><td>60.1</td><td>66.4</td><td>8.51</td><td>0.11</td></tr><tr><td>Fe</td><td>41.1</td><td>45.2</td><td>36.9</td><td>3.29</td><td>0.19</td></tr><tr><td>Mn</td><td>2.64</td><td>2.26</td><td>2.51</td><td>0.10</td><td>0.14</td></tr><tr><td>Mo</td><td>1.16</td><td>1.06</td><td>1.11</td><td>0.07</td><td>0.65</td></tr><tr><td>Se</td><td>0.48</td><td>0.42</td><td>0.50</td><td>0.02</td><td>0.16</td></tr><tr><td>Sn</td><td>0.73</td><td>0.78</td><td>0.66</td><td>0.05</td><td>0.33</td></tr><tr><td>Zn</td><td>33.2</td><td>25.6</td><td>27.9</td><td>2.11</td><td>0.12</td></tr><tr><td colspan="6"><strong>Macro-element, mg/kg</strong></td></tr><tr><td>Mg</td><td>131</td><td>128</td><td>131</td><td>3.8</td><td>0.59</td></tr><tr><td>Ca</td><td>40</td><td>41</td><td>37</td><td>1.7</td><td>0.30</td></tr><tr><td>K</td><td>2450</td><td>2482</td><td>2533</td><td>74.6</td><td>0.70</td></tr><tr><td>Na</td><td>684</td><td>707</td><td>721</td><td>18.7</td><td>0.46</td></tr><tr><td>P</td><td>2655</td><td>2693</td><td>2692</td><td>72.0</td><td>0.74</td></tr><tr><td>S</td><td>1555</td><td>1506</td><td>1561</td><td>31.9</td><td>0.33</td></tr></tbody></table>

SEM: Mean standard error of the means.

doi:10.1371/journal.pone.0167747.t003

respective treatments (S2 Fig). However, interclass PCA did not identify any gene having extreme count profiles, which may have contributed to the transcriptome dispersion of animal BC7_Conv with respect to its treatment group. This sample was therefore retained for all the analyses. The total amount of expressed genes in liver was similar in all dietary treatments (Conv = 15,503, LowI = 14,994; Nutr = 15,576 genes). A total of 14,549 common genes were expressed in all groups. The correlation of mean gene expression levels between groups was very high (r = 0.92), suggesting that the major fraction of the liver transcriptome was not affected by nutrient ingestion.

**Analysis of alternative splicing.** We used the RNA-sequencing data to characterize and compare the patterns of alternative splicing in the 3 dietary groups. Multivariate Analysis of Transcript Splicing (rMATS) identified a total of 266 differential alternative splicing events: 14 mutually exclusive exons (MXE), 205 skipped exons (SE), 14 alternative 5' splice sites (A5SS) and 32 retention introns (RI), while alternative 3' splice sites (A3SS) events were absent. All 266 AS events were distributed across 238 multi-exon genes. When testing for differential alternative splicing in all comparisons, a total of 9 significant splicing events in 7 genes were found after multiple testing correction (Table 4, S3 Table). Interestingly, 2 SE in fibronectin 1 (FN1) (exon 33 and exon 25), 1 SE in leucine rich repeat (In FLII), interacting protein 1 (LRRFIP1) (exon 14), and 1 RI in complement component 2 (C2) were significantly associated with the Nutr diet in comparison with both LowI and Conv (Table 4).

**Identification of differentially expressed gene.** RNA-sequencing data were used to study the effect of dietary composition on the liver transcriptome, with a focus on the identification of genes whose expression was significantly altered by LowI relative to Conv, Nutr relative to LowI, and Nutr relative to Conv. We identified 89 differentially expressed genes (DEGs), of which 52 were over-expressed and 37 were under-expressed in LowI compared with Conv.

<footnote n="1">¹Mean values are reported for each diet.</footnote>

<footnote n="2">²P-values after Kruskal Wallis non-parametric test.</footnote>
MDS 1x2

<figure><img src="image_5.png" /><figcaption>Fig 1. Multidimensional scaling (MDS) analysis using log2 of all normalized unique gene counts (n = 24,596). Conv: conventional diet; Lowl: low-impact diet; Nutr: nutraceutical diet.</figcaption></figure>

doi:10.1371/journal.pone.0167747.g001

<table><caption>Table 4. Significant differences in splicing events among the three diet groups (conventional, low-impact and nutraceutical).</caption><thead><tr><th>GeneID</th><th>GeneSymbol</th><th>Comparison</th><th>Event Type</th><th>FDR</th><th>IncLevelDifference</th></tr></thead><tbody><tr><td>ENSBTAG00000011424</td><td>TPM2</td><td>Lowl vs Conv</td><td>MXE</td><td>0.0152</td><td>0.045</td></tr><tr><td rowspan="4">ENSBTAG00000008300</td><td rowspan="4">FN1</td><td>Nutr vs Conv</td><td>SE</td><td>0.0006</td><td>0.019</td></tr><tr><td>Nutr vs Conv</td><td>SE</td><td>0.0204</td><td>0.024</td></tr><tr><td>Nutr vs Lowl</td><td>SE</td><td>0.0162</td><td>0.017</td></tr><tr><td>Nutr vs Lowl</td><td>SE</td><td>0.0162</td><td>0.022</td></tr><tr><td rowspan="2">ENSBTAG00000005354</td><td rowspan="2">LRRFIP1</td><td>Nutr vs Conv</td><td>SE</td><td>0.0031</td><td>0.113</td></tr><tr><td>Nutr vs Lowl</td><td>SE</td><td>0.0162</td><td>0.124</td></tr><tr><td rowspan="2">ENSBTAG00000007450</td><td rowspan="2">C2</td><td>Nutr vs Conv</td><td>RI</td><td>0.0272</td><td>0.030</td></tr><tr><td>Nutr vs Lowl</td><td>RI</td><td>0.0001</td><td>0.033</td></tr><tr><td>ENSBTAG00000006835</td><td>MCAM</td><td>Nutr vs Lowl</td><td>SE</td><td>0.0006</td><td>0.284</td></tr><tr><td>ENSBTAG00000021319</td><td>MFF</td><td>Nutr vs Lowl</td><td>SE</td><td>0.0180</td><td>-0.061</td></tr><tr><td>ENSBTAG00000014791</td><td>CTH</td><td>Nutr vs Lowl</td><td>RI</td><td>0.0139</td><td>-0.202</td></tr></tbody></table>

Lowl: low-impact diet; Conv: conventional diet; Nutr: nutraceutical diet

FDR: false discovery rate

MXE: mutually exclusive exon; SE: skipped exon; RI: retained intron; IncLevelDifference: difference in the exon inclusion level (%) between the experimental groups tested in each comparison

doi:10.1371/journal.pone.0167747.t004
Comparing Nutr and LowI, we identified a total of 166 DEGs (79 up-regulated and 87 down-regulated). Regarding Nutr relative to Conv, 113 genes were differentially expressed (57 up-regulated and 56 down-regulated). Only three DEGs were commonly regulated among all 3 groups: nocturnin (CCRN4L or NOCT), leucine rich adaptor protein 1-like (LURAP1L), and pyruvate dehydrogenase kinase isozyme 4 (PDK4) (S3 Table).

**Gene ontology and pathway analysis.** We performed gene set enrichment analysis in order to identify over-represented gene ontology (GO) terms and Kyoto encyclopedia of genes and genomes (KEGG) pathways. In all comparisons, DEGs were strongly associated with lipid metabolism, molecular transport, and small molecule biochemistry. The full list of significantly enriched pathways is given in S4 Table. In addition, we characterized the transcription factors and microRNAs (miRNAs) that regulate DEGs (S5 Table).

The 89 DEGs found when comparing LowI with Conv were strongly associated with the isoprenoid biosynthetic process and glycogen metabolic process (S3 Table). In addition, the cholesterol biosynthesis canonical pathways and LXR/RXR activation were overrepresented (S4 Table). Upstream regulator analysis by Ingenuity Pathway Analysis (IPA) showed that the RAR-related orphan receptor C (RORC) gene, together with the insulin induced gene 1 (INSIG1), sterol regulatory element binding transcription factor 1 (SREBF1), sterol regulatory element binding transcription factor 2 (SREBF2), and peroxisome proliferator-activated receptors alfa (PPARA) played a central role in mediating the cholesterol biosynthesis in animals fed with the LowI diet (S5 Table). In particular, these TFs affected the expression of genes involved in cholesterol biosynthesis (S5 Table). The 166 DEGs found when comparing Nutr with LowI treatments were related to metabolism of glycogen and carbohydrates, lipid homeostasis together with cholesterol biosynthesis and triglyceride metabolism, as well as fat cell differentiation and steroid biosynthetic process (S3 Table). The list of top upstream regulators revealed by IPA included INSIG1 followed by the cooperatively transcriptional cofactors RORC, insulin receptor substrate (IRS) 2 and FA synthase (FASN) (S5 Table). TFs act in regulatory networks and can drive the expression of other TFs in a feed-forward and feedback manner. In fact, these TFs were in turn regulated by SREBF1 and SREBF2 together with other genes markedly down-regulated and mainly involved in cholesterol biosynthesis (S5 Table).

When evaluating the effects of Nutr relative to Conv diet (113 DEGs), we observed that the most significantly enriched GO terms regarded the metabolism of energy reserve and differentiation of fat cell (S3 Table), including gluconeogenesis, triglyceride metabolic process, cholesterol biosynthesis and lipid homeostasis. Enrichment of the glutathione metabolism pathway was also found, as was activation of TR/RXR (S4 Table). Top upstream regulators included genes associated with the immune system, such as tumour necrosis factor (TNF) and the interleukin 1B (IL1B) (S5 Table; Fig 2), as well as TF playing important roles in lipid metabolism (i.e. the down-regulated INSIG1, IRS2 and the up-regulated and RORC (Fig 2). In addition, Nutr down-regulated the expression of several target genes for PPARA compared to Conv diet (S5 Table).

## Correlations between differentially expressed genes and phenotypic traits

Pearson correlations between DEGs and phenotypic traits are reported in Fig 3A–3C. Moderate to strong significant correlations ($0.62 < r < 0.73$) were found between DEGs and *in vivo* fatness scores (Fig 3A). Most of the DEGs positively correlated with the *in vivo* fatness scores were involved in the regulation of energy metabolic process (IRS2), regulation of cholesterol biosynthesis (INSIG1), and circadian clock (e.g. NOCT and NR1D1). In contrast, the RORC and glycerol-3-phosphate acyltransferase mitochondrial (GPAM) genes, as well as the
<figure><img src="image_6.png" /><figcaption>**Fig 2. Regulatory network of transcription factors, miRNAs and target genes using Cytoscape.** Schematic visualization of shared associations between transcription factors (TFs), miRNAs, and differentially expressed genes when comparing Nutr with Conv. The network is displayed graphically as nodes (genes, TFs, and miRNAs) and edges (the biological relationship between nodes). Node color intensity indicates the expression of the association: red = upregulation with Nutr diet, green = downregulation with Nutr diet. Node shape indicates whether it is a TF (triangle) or an miRNA (round) or other kinds of molecule.</figcaption></figure>

doi:10.1371/journal.pone.0167747.g002

glycoprotein (transmembrane) Nmb (GPNMB), which is associated with long-chain FA (LCFA) transport, were negatively correlated with fatness score ($r = -0.732, -0.90$ and $-0.74$, respectively). Additionally, a large number of DEGs correlated with pH in liver (Fig 3B). For instance, genes involved in the regulation of cation transmembrane transport [e.g. fibroblast growth factor 12 (FGF12) and kell blood group metallo-endopeptidase (KEL)] showed moderate positive correlations with liver pH ($r = 0.62$ and $0.71$, respectively).

The association between liver FA composition and gene expression was calculated as a means of potentially identifying the association between the expression of lipid metabolism-related genes and FA contents in liver. Using the Pearson correlation approach it was possible to split liver FAs into two main clusters, which differed largely in their overall gene expression profiles. Gene cluster 1 contained 77 genes and had moderate to strong negative correlations
<figure><img src="image_7.png" /><figcaption>Fig 3. Heatmap of correlations between gene expression changes and phenotypic variations. Correlations were estimated for growth performances and carcass traits (A), liver weight, pH, and chemical composition (B), and liver fatty acid composition (groups of fatty acids) (C).</figcaption></figure>

The heatmap was generated using the CorrPlot package in R. The gene was included in the heatmap when a significant correlation for at least one phenotype was found (*P* < 0.05, Pearson product-moment correlations). Blue indicates negative correlations with values *r* < -0.51; grey indicates correlations values in the range -0.50 *r* < 0.50; red color indicates positive correlations with *r* values > 0.51. BW: body weight; DG: daily gain; FS_IV: fat covering *in vivo*; MS_IV: body conformation; CW: carcass weight; MS_PM: carcass muscularity score; FS_PM: carcass fatness score; WT: liver weight. DM: dry matter; SFA: saturated fatty acids; MUFA: monounsaturated fatty acids; PUFA: polyunsaturated fatty acids; SCFA: fatty acids < C16; C16: fatty acids with 16 carbon chains; LCFA: fatty acids > C16; OCFA: odd-chain fatty acids; BCFA: branched-chain fatty acids; TFAC18_mono: Trans C18:1 fatty acids; CLAct_tc: conjugated linoleic acid *c*, *t*, *t*, *c* isomers; n_Ratio: ratio n-6/n-3.

doi:10.1371/journal.pone.0167747.g003

with n-6 FA, n-6/n-3 ratio, FAs with 16 carbon chains, MUFAs, and SFAs, moderate to strong positive correlations with n-3 FA, CLA *c*, *t*, *t*, *c*, *trans* 18:1, as well as PUFAs and LCFAs, and no specific correlation with odd-chain FAs, branched-chain FAs and short-chain FAs (SCFAs) (Fig 3C). Gene cluster 2 contained 68 genes and had the opposite correlations. That is, we found CLA and n-3 to have moderate to strong correlations with genes related to the regulation of cholesterol biosynthesis (*GPAM*, *INSIG1*), glycogen biosynthetic process (protein phosphatase 1 regulatory (*PPP1R*) 3B, and *PPP1R3C*), circadian clock [e.g. aryl hydrocarbon receptor nuclear translocator-like (*ARNTL*), *NOCT*, and *NR1D1*)], and regulation of FA metabolism (*INSIG1*, *IRS2*, and *PDK4*). On the other hand, SCFAs and FAs with 16 carbon chains showed moderate positive correlations with *INSIG1* (*r* = 0.71), *IRS2* (*r* = 0.54), and *FASN* (*r* = 0.67). In addition, n-6 FA was found to correlate with genes associated with carbohydrate catabolic process (e.g. glucose-6-phosphatase, catalytic subunit (*G6PC*), *r* = 0.6; myoinositol oxygenase (*MIOX*), *r* = -0.77), and the sphingolipid-mediated signaling pathway (sphingosine-1-phosphate receptor 3 (*S1PR3*), *r* = 0.67; spinster homolog 2 (*SPNS2*), *r* = -0.69). Interestingly, CLA, n-3, and n-6 FA also presented moderate to strong correlations with genes involved in inflammatory response, such as *LURAP1L* and insulin-like growth factor binding protein 1 (*IGFBP1*), and with stress-response genes, such as *DDIT4* and *BTG* family member 2 (*BTG2*). Total PUFAs were highly associated with the expression of genes involved in LCFA transport (*GPNMB*; *IRS2* and peripilin 2 (*PLIN2*)).
Regarding the correlation of mineral contents with the expression profile of the DEGs, the highest number of significant and relevant correlations were found with Ca, Fe, Mn, Se and Zn (S6 Table). Significant and moderate to strong correlations were found for some minerals and genes involved in the cellular transport system: e.g., Mn and Zn were negatively correlated with ATP-Binding Cassette, sub-family C member 5 (ABCC5) (r = -0.70 and -0.83, respectively), and Ca and Se were negatively correlated with ATP-Binding Cassette, Sub-Family G Member 8 (ABCG8) (r = -0.62 and -0.77, respectively). Positive moderate correlations were found for Zn contents and metallothionein (MT) 1A (r = 0.64), 1E1 (r = 0.78) and 1E2 (r = 0.68).

## Discussion

In the present experiment, we evaluated the effects of feeding strategies designed to have the least interference with growth performances but differing in the supply of specific nutrients in young Belgian Blue × Holstein bulls. As expected, the dietary changes had small or no relevant influence on growth performances and main carcass traits, and on the composition of major nutrients in the liver, but they had a major influence on liver FA profile and clearly affected the transcriptome expression of bovine liver.

## Dietary treatments and phenotypic traits

The use of the LowI diet, where extruded full-fat soybeans replaced the soybean meal of the Conv diet, notably decreased protein and increased n-6 FA intakes. Inadequate protein intake may impair growth and host immunity with particularly detrimental effects on the T cell system, resulting in an increased incidence of infections [22], and may also affect the insulin/glucagon ratio [23] and the rate of hepatic glycogen synthesis [24]. However, an excess of protein can also be detrimental to ruminant health, as an excess of rumen-degradable N increases lameness and laminitis, and certain types of proteins may cause allergic histaminic reactions [25]. The effect of including extruded full-fat soybean rich in 18:2c9,c12 in the diet was evident primarily in an increase in the proportion of 18:2c9,c12 in the liver fat. Other experiments conducted on a larger number of bulls found that increased amounts of extruded full-fat soybeans increased the 18:2c9,c12, 18:2c9,t11, 18:1t11, or other intermediates of rumen biohydrogenation, in intramuscular, intermuscular and subcutaneous fat in beef [25-26].

Notably, Nutr increased n-3 FAs, CLA, and vitamin E intakes compared with LowI, due to the inclusion of linseed CLA, and vitamin E. This feeding treatment had no effect on growth performances, except for a reduction in the *in vivo* subcutaneous fat covering. Data on the effects of CLA on body fat conflict, as the reduction of body fat was not always confirmed in growing ruminants [25-26, 27]. Notable alteration to the FA profile of liver fat was instead induced without apparently affecting the liver fat content. In particular, linseed, rich in 18:3c9, c12,c15, caused a marked increase in this FA in the liver fat. In general, the liver fat of bulls fed Nutr had higher proportions of LCFA, n-3 FA, and CLA isomers (2.5-fold) compared with LowI. Compared with Conv, Nutr and LowI caused a significant increase in the liver content of rumen hydrogenation intermediates, such as *trans*-18:1 FA and CLA, which were found to inhibit *de novo* FA synthesis in the mammary gland [28].

Concerning mineral concentrations in liver, we observed that most of the trace and macro-mineral contents of cattle liver in our study were in line with data from the literature [29-31]. Changes in the dietary content of macro-elements in the LowI and Nutr diets only nominally affected the liver contents compared with Conv. A larger number of animals would be required to draw statistical conclusions from these differences.
Effects of reduced protein and mineral supplementation, and increased
soybean oil supplementation on bovine liver transcriptome

The changes in the protein and lipid profiles characterizing LowI diet were associated with down-regulation of glycogenesis and up-regulation of cholesterol biosynthesis in bovine liver. Liver metabolism can be regulated by various metabolites, e.g. NEFAs and glucose [32]. The LowI diet seemed to affect the regulation of insulin homeostasis, as suggested by the over-expression of insulin target genes, such as FASN and PDK4 (S5 Table). Coupled to that, up-reg-ulation of target genes for SREBF1, SREBF2 and PPARA supported the hypothesis that LowI enhance regulation of cholesterol biosynthesis in bovine liver. The activation of RORC could also play a role on the induction of cholesterol biosynthesis, the inhibition of glycogen synthe-sis and the regulation of the circadian clock. Circadian clock controls energy metabolism through several genes involved in FA, cholesterol and glucose metabolism in liver, which are regulated in circadian or diurnal patterns [33]. However, our results did not allow ascertaining if the ingestion of higher amounts of PUFA caused the increase in cholesterol biosynthesis and consequently, an enhanced biliary cholesterol secretion [34]. Besides to the aforementioned changes in energy metabolism, we observed the activation of genes involved in immune response which may be explained by the lower intake of protein and/or the higher intake of dietary n-6 PUFAs (linoleic or arachidonic acid) which alter cell membrane fluidity and indi-rectly affect inflammatory processes in that they increase the pro-inflammatory response [35].

Over the past decade, miRNAs have emerged as novel elements in the rapid, reversible reg-
ulation of transcription and translation [36]. MiRNA are small, non-coding RNA molecules
(~19–24 bp in length) that are synthesized from short hairpin precursors and that reportedly
degrade or inhibit translation of their target genes by binding to the 3' untranslated region
(UTR) of coding mRNA [36]. Interestingly, the LowI diet modified the expression of two mir-
185 target genes, which are involved in cholesterol biosynthesis (S5 Table), suggesting that this
miRNA may also play a role in liver lipid metabolism along with the above-mentioned TFs.
Indeed, miR-185 is a key micro-RNA identified as a regulator of de novo cholesterol biosynthe-
sis and low-density lipoprotein uptake; inhibition of miR-185 resulted in decreased SREBP-2-
dependent gene expression, LDL uptake, and HMG-CoA reductase activity [37].

Effects of n-3, CLA and vitamin E supplementation on bovine liver transcriptome

A large degree of transcriptomic adaptation to the n-3 and CLA FAs was observed in the liver of bulls fed the Nutr diet. Biological function analysis revealed under-expression of TFs regulating lipid metabolism, such as *INSIG1* and *IRS2* (Fig 2). In addition, FASN and the target genes for *SREBF1* and *SREBF2* were strongly inhibited by the Nutr diet (Fig 2). When cells have sufficient sterol levels, *INSIG1* retains the *SREBP1* cleavage-activating protein (SCAP)-*SREBP1* in the endoplasmic reticulum and consequently inhibits *SREBP1*-mediated gene expression [38]. In addition, *INSIG1* is itself a target of *SREBP*, and is repressed when *SREBPs* are inactive [39]. In bovine mammary epithelial cells, down-regulation of *INSIG1* and *FASN*, a key enzyme involved in *de novo* FA synthesis by oleic acid, linoleic acid, VA and alpha-linoleic acid (n-3), was observed [40–41]. Hiller et al. [42] also found evidence that supplementation with n-3 FA caused reduced expression of various genes involved in FA synthesis, including *FASN*, in Holstein cow liver. In light of the above observations, it is tempting to speculate that PUFAs (in particular CLA and n-3 FA) reaching the liver affect the expression of both *SREBP1* and *INSIG1* to inhibit *SREBP1*-mediated gene expression, and consequently, at least partially, likely reduce the rate of lipogenesis occurring in bovine liver in a negative feedback mechanism. In addition, the lower expression of *IRS2* may support an increase in the
extent of FA oxidation [43] in Nutr, compared with other diets. With Nutr, the increased amount of PUFAs reaching the liver likely increased hepatic oxidation of FA, driving liver metabolism towards lipid uptake and utilization, unlike what was observed in bulls fed LowI. This hypothesis is further supported by the down-regulation of *PDK4* (*S3 Table*). Depletion of *PDK4* increases PDC activity, which allows pyruvate to be channeled to the tricarboxylic acid cycle for complete oxidation and as a result pyruvate is not available for gluconeogenesis [44]. It would seem that under-expression of *PDK4* came about by an increased amount of PUFAs reaching the liver after administration of Nutr, consistently with previous studies [45]. Additionally, down-regulation of the key gluconeogenic factor G6PC seems to be consistent with a negative effect on gluconeogenesis in liver. Considering the effects of Nutr on genes related to immune/stress response, our results seemed to support the hypothesis that supplementation with n-3 FA and CLA may decrease the production of inflammatory eicosanoids, cytokines, and ROS, and also decrease *de novo* FA synthesis in liver, as published elsewhere [46]. This could explain why CLA supplementation may improve the health status of ruminants [21].

Several target genes for the miR Let-7 and Let-7a-5p were down-regulated in the animals on the Nutr diet (*S5 Table*), supporting previous evidence concerning the role of the Let-7 family of miRNA in regulating glucose metabolism in the liver [47]. Target genes for miR-21 were differentially expressed in Nutr vs LowI (6), and Nutr vs Conv (3) (*S5 Table*), according to recent findings concerning its putative effect on the regulation of triglyceride and cholesterol metabolism [48].

The role of dietary nutrients, such as carbohydrates, FAs, and protein, in regulating alternative splicing events and in controlling nutrient metabolism through variations in the splicing of the pre-mRNA encoding enzymes involved in various metabolic pathways has been recently highlighted. For instance, splicing of the pre-mRNA encoding *G6PD* is modulated by carbohydrates and FA [49]. Interestingly, we found an SE event in *PFKFB2*, a gene involved in the control of glycolysis, which was also upregulated in the Nutr group (*S3 Table*). However, the differences in the inclusion levels of this splicing event between our experimental groups were not significant (data not shown). On the other hand, in the Nutr group we found a significant difference in the percentage of inclusion of AS events in genes involved in the innate immune response, such as *FN1*, *LRRFIP1* and *C2* [50–52]. The possibility that these events may be induced in the individuals on the Nutr diet, favoring alternatively spliced isoforms with specific roles, cannot be excluded. Interestingly, in recent years multiple roles of alternatively spliced FN isoforms have been identified in humans [53].

## Correlations between gene expression changes and phenotypic variations

The highest number of significant and relevant correlations with DEGs was found with those traits that were significantly different in at least one comparison between diets (Fig 3A-3C). Interestingly, fat covering *in vivo* showed a moderate positive correlation with *LURAP1L* (*r* = 0.79), which was one of the three genes differentially expressed in all comparisons. Its paralog, *LURAP1*, was shown to play a crucial role in obesity-induced inflammation [54]. As expected, moderate to strong correlations were found between fat covering *in vivo* and genes related to energy reserve metabolic process. A moderate positive correlation was also found between fat covering *in vivo* and *NR1D1* (*r* = 0.73), which has a role in the network of Fig 2, further confirming the link between energy metabolism and the expression of clock genes [55].

Genes involved in membrane transport showed moderate positive correlations with liver pH, providing support for a connection between hepatic membrane transport and pH regulation mechanisms [56].
Finally, the results of correlation analyses of liver FA profiles revealed that CLA and n-3 FA might exert similar effects on bovine liver fat metabolism (increased FA oxidation and decreased lipid synthesis), as they seemed to act on lipid-related genes (e.g., key TFs, such as *INSIG1* and *IRS2* together with some of their target genes) with likely equivalent effects (same degree and direction of correlations) (Fig 3C). On the other hand, it is interesting to note that, when looking at the correlation values of the DEGs involved in the inflammatory/stress response, such as *LURAP1L*, *BTG2*, and *IGFBP1*, CLA and n-3 FA seemed to have the opposite effect to that of n-6 FA (negative vs positive correlations), thus confirming the anti-inflammatory properties of the former and the pro-inflammatory potential of dietary n-6 PUFA [46].

In general, weak correlations between DEG and mineral contents were found (S6 Table), suggesting that the variation in the dietary mineral content probably did not have a great effect on the bovine liver gene expression profile. Interestingly, moderate positive correlations were found between *MT1A* and *MTE1* and liver Zn content, confirming that cells are likely to protect themselves from zinc toxicity by inducing MTs, which bind it tightly by sequestering it in organelles, or by exporting it [57].

## Conclusions

In summary, this trial suggested the feasibility of using low environmental impact feeding regimes (low protein and mineral contents) in livestock production without compromising animal performances and carcass characteristics. The RNA-sequencing data here obtained have broadened our understanding of the effects of variations in feeding strategies on bovine liver transcriptome, and have unraveled, at least partially, the mechanisms underlying the relationships between nutrients and phenotypic variations. In particular, the balance between glucose and lipid metabolism was affected by both low-impact diets (LowI and Nutr), while the greatest effects on liver transcriptome resulted from the addition of linseed and CLA (Nutr). *INSIG1*, *IRS2* and *RORC*, key TFs controlling the expression of lipid-related genes, clearly played a pivotal role in mediating the observed effects on gene expression. On the other hand, the changes in mineral composition seemed not to considerably affect bovine liver gene expression profile. The increased supplementation with vitE occurring in Nutr could likely be related to the regulation of glutathione metabolism.

Thereby, our results helped in clarifying the role of bovine liver in the biological response to specific nutrients variations. The regulation of liver metabolic processes involved in lipid metabolism may subsequently influence fatty acid deposition into the fat of ruminant products. It is clear however, that validation on a larger set of biological samples is needed in order to confirm the present RNA-sequencing findings and draw more robust conclusions about the interactions between diet and bovine liver metabolism.

## Material and Methods

### Animals, housing and diets

This study was part of a large experiment aimed at investigating the effects of different feeding strategies on animal performances, carcass traits, and chemical and FA compositions of different tissues (muscle, liver, and adipose tissue) in 48 crossbred heifers and young bulls (15–19 months old at slaughter). The project was approved by the Ethical Committee for the Care and Use of Experimental Animals of the University of Padua (CEASA) and all methods were carried out in accordance with the approved guidelines. Following the experimental protocol for animal care, health status was monitored daily by a technician and three times per week by a veterinarian.
<table><caption>Table 5. Ingredient and chemical composition of the diets.</caption><thead><tr><th></th><th>Conventional diet (Conv)</th><th>Low impact diet (Lowl)</th><th>Nutraceutical diet<fnref n="1" /> (Nutr)</th></tr></thead><tbody><tr><td><strong>Ingredient composition, kg/d:</strong></td><td></td><td></td><td></td></tr><tr><td>Corn silage</td><td>2.53</td><td>2.53</td><td>2.53</td></tr><tr><td>Corn meal</td><td>3.61</td><td>3.61</td><td>3.61</td></tr><tr><td>Wheat bran</td><td>0.62</td><td>0.62</td><td>0.62</td></tr><tr><td>Wheat straw</td><td>0.65</td><td>0.65</td><td>0.65</td></tr><tr><td>Sugar beet pulp</td><td>1.06</td><td>1.06</td><td>1.06</td></tr><tr><td>Grape-seed meal</td><td>0.18</td><td>0.18</td><td>0.18</td></tr><tr><td>Soybean meal (44% CP)</td><td>1.29</td><td>-</td><td>0.20</td></tr><tr><td>Corn meal</td><td>-</td><td>0.79</td><td>0.53</td></tr><tr><td>Extruded soybeans, full fat</td><td>-</td><td>0.48</td><td>-</td></tr><tr><td>Linseed, full fat</td><td>-</td><td>-</td><td>0.48</td></tr><tr><td>Vitamins-minerals</td><td>0.20</td><td>0.20</td><td>0.20</td></tr><tr><td><strong>Chemical composition:</strong></td><td></td><td></td><td></td></tr><tr><td>Dry matter, g/kg as fed</td><td>757</td><td>755</td><td>757</td></tr><tr><td>Crude protein, g/kg DM</td><td>144</td><td>107</td><td>107</td></tr><tr><td>Neutral detergent fiber, g/kg DM</td><td>300</td><td>295</td><td>296</td></tr><tr><td>Starch, g/kg DM</td><td>360</td><td>418</td><td>401</td></tr><tr><td>Ether extract, g/kg DM</td><td>30</td><td>42</td><td>55</td></tr><tr><td>Ash, g/kg DM</td><td>44</td><td>39</td><td>41</td></tr><tr><td>Vitamin A</td><td>4123</td><td>4207</td><td>4176</td></tr><tr><td>Vitamin E</td><td>25</td><td>30</td><td>116</td></tr><tr><td><strong>Net energy content, UFC/kg DM</strong></td><td>1.02</td><td>1.03</td><td>1.04</td></tr></tbody></table>

doi:10.1371/journal.pone.0167747.t005

Briefly, 24 Belgian Blue × Brown Swiss (12 ♂ and 12 ♀), and 24 Belgian Blue × Holstein (12 ♂ and 12 ♀) animals were weighed and housed in 12 pens (4 animals per pen), balanced by sex, breed, age and initial BW, and considering the representativeness of half-sib families based on pedigree information. The animals were assigned to 1 of 3 dietary treatments following a facto-rial arrangement. All animals belonging to the same pen were fed the same dietary treatment.

After an adaptation period, a first group (Conv) was fed a conventional diet containing 14.4% crude protein (CP) on a DM basis, mainly from soybean meal, and supplementation of minerals and vitamins for the whole experimental period of 9 months (Table 5). The second group was fed for the whole experiment a LowI diet containing 10.7% CP on a DM basis, obtained by substituting all the soybean meal with a small amount of full-fat soybean, and reducing the mineral supplementation (Table 5).

In the third treatment, the animals were fed the LowI diet for the first 5 months, and the Nutr diet for the following 4 months (for young bulls). In the Nutr diet, the full-fat soybeans in the LowI diet were replaced with linseed rich in n-3 FA, adding a high level of vitamin E [as an antioxidant to compensate for a possible reduction in beef shelf-life due to increased oxidation of the PUFAs derived from linseed], and top dressed with 80 g/d of a commercial CLA mixture (SILA, Noale, Italy). Whole linseed was chosen because this oilseed is rich in n-3 FA (54.2% alpha-linolenic acid), and its seed coat may protect PUFA from rumen biohydrogenation and increase the passage of PUFA to the duodenum. The diets were distributed once a day at 8.00 a.m., and the animals had free access to the manger and drinkers. The amount of each feed

<footnote n="1">¹The Nutraceutical diet was also top dressed with 80 g/d of a commercial rumen protected conjugated linoleic acid product (SILA, Noale, VE) providing 5.57 g/d of *cis-9*, *trans-11* CLA and 5.41 g/d of *trans-10*, *cis-12* CLA. All diets contained 0.79 g/kg dry matter of Na, and 401 IU/kg dry matter of Vitamin D3. CP: crude protein; DM: dry matter.</footnote>
ingredient loaded into the mixer-wagon and the weight of the mixture uploaded into the man-
ger of each pen were recorded daily. The orts remaining in the mangers were weighed and
sampled by pen weekly.

Animal BW and ADG were recorded for each bull 3 times during the experiment. An oper-
ator licensed for carcass evaluation according to the SEUROP grading system (European
Commission, 1991) made monthly evaluation of the *in vivo* body condition of each bull as
described by Schiavon et al. [21]. Conformation was linearly scored from S+ (all muscle pro-
files extremely convex; superior muscle development) to P- (all muscle profiles concave to
very concave; very poor muscle development). Conformation was expressed in numerical
terms: S+ = 6.33; S = 6; S- = 5.66; . . . , P+ = 1.33; P = 1.00; P- = 0.66. Fat covering was scored
linearly from 1 (very lean) to 5 (very fat) using a combined visual and palpation approach, tak-
ing into account the presence and thickness of subcutaneous fat depots at specific points of the
body.

Slaughtering, carcass traits, and sample collection

Animals were transported to the slaughterhouse after 280 days on experimental feed with an average target slaughtering weight of approximately 670 kg. They were all fasted for 1 day and slaughtered all together at 532 ± 36 days old. The animals were killed using captive bolt pistol in the stunning process and all efforts were made to minimize suffering. Carcasses were indi- vidually weighed and scored for conformation and fat covering according to the SEUROP sys- tem (European Commission, 1991). Dressing percentage was computed as the ratio between the carcass weight 24 h after slaughter and BW.

Immediately after slaughter, liver tissue samples were collected from all the animals,
weighed, then stored at -20° for chemical and FA analyses. The liver tissue samples for RNA-
sequencing analyses (~30 mg) were placed in cryogenic, sterile, DNAase- and RNAase-free
vials (Cryo.s, Greiner Bio One, Frickenhausen, Germany) and immediately snap-frozen in liq-
uid nitrogen and stored at -80°C until RNA extraction. Since sex and breed differences in liver
transcriptome have been reported for several species, including cattle [58], only the 12 males
from the same crossbreed (Belgian Blue × Holstein) were selected for subsequent analyses.
Female animals were excluded because of the lower importance of heifers for meat production
and as hormonal fluctuations might affect liver gene expression profiles in mammals [59].
Belgian Blue × Holstein crossbreed young bulls were instead chosen respect to Belgian
Blue × Brown Swiss crossbreds because Holstein breed is much more used than Brown Swiss
in dairy farms worldwide and as they were more balanced in terms of half-sib families. Data
for the transcriptome profiling therefore came from 4 Belgian Blue × Holstein bulls for the
Conv diet, 4 for LowI, and 4 for Nutr. One bull on the Conv diet was removed from the study
due to health problems unrelated to the dietary treatment.

Chemical analysis

**Liver proximate composition.** The 50 g liver samples were sliced, ground and homoge-
nized using a food processor (Grindomix GM200; Retsch, Haan, Düsseldorf, Germany), and
subsamples of 20 g were used to determine their chemical composition according to AOAC
(2000). In brief, moisture was determined after drying at 102°C for 16 h, while ash was deter-
mined after mineralization at 525°C. Total lipids were analyzed by extraction with petroleum
ether, and protein was estimated by difference, in accordance with Schiavon and colleagues
[25]. Liver pH was measured on fresh samples using a Delta Ohm HI-8314 pH-meter (Delta
Ohm, Padova, Italy).
**Liver fatty acid profile.** The left lobe of each liver was sliced, ground, mixed, and then homogenized for 10 sec at 4500 g (Grindomix GM200; Retsch, Haan, Düsseldorf, Germany). Two subsamples of 50 to 60 g each were used to determine the FA composition. Samples were stored at -20°C until analysis. After conservation, a subsample (2.00 ± 0.005 g of liver, approximately 44 mg of fat) was freeze-dried with a CoolSafe 90-80 freeze dryer (Scanvac, Stockholm, Sweden), according to the manufacturer's instructions. Freeze-dried samples were methylated as described by Jenkins [60] to limit trans-isomerization of the FAs [61]. The samples obtained were analyzed for their FA profiles using a GC × GC instrument (7890A; Agilent Technologies, Santa Clara, CA, USA) with two columns in series and equipped with a modulator (G3486A CFT; Agilent Technologies), an automatic sampler (7693A; Agilent Technologies) and a flame ionization detector (FID) connected to the Agilent Chemstation chromatography software (Agilent Technologies). Operating conditions and GC reference standards were previously reported in detail [27]. The liver FA composition was expressed as grams per 100 g of total FA.

**Liver mineral profile.** Liver samples were digested using the microwave digestion method (Milestone Start D, Milestone Srl, Sorisole, BG, Italy) (AOAC, 2000). Samples of approximately 1.0 g were digested with 7 ml of HNO₃ and 2 ml of H₂O₂ 30% in the microwave digestion system. The samples and acid mixture were placed in suitable inert polymeric microwave vessels. Each vessel was sealed and heated in the microwave digestion system. The temperature program was as follows: 15 min to reach 200°C and held for 18 min. The resulting solution was cooled to 35°C and diluted to 25 ml with distilled water. Determination of 27 micro- and macro-elements in this clear solution was carried out by SPECTRO ARCOS ICP-OES (SPECTRO Analytical Instruments, GmbH, Kleve, Germany). Instrument operating parameters were optimized for aqueous solutions with undissolved organic material. Calibration standards were matched with 1% absolute ethanol (Prolabo VWR International PBI srl, Milan, Italy). The concentration ranges of the calibration solutions were 0-250 mg/L for macro-elements and 0-20 mg/L for micro-elements. The elements to be determined were added from single element solutions (Inorganic Ventures, Christiansburg, VA, USA). The measured and certified values of all elements were in excellent agreement. Element contents were expressed as mg/kg fresh tissue.

## Statistical analyses of phenotypic data

A Kruskal Wallis non-parametric test with Bonferroni correction was used to identify differences in the phenotypic traits among the 3 dietary treatments (Conv, LowI, Nutr). The level of significance was set at an adjusted P-value < 0.05.

## Liver transcriptome profiling

**RNA extraction.** Total RNA was isolated from liver samples (30 mg) stored at -80°C using the RNAeasy Mini kit (Qiagen, Hilden, Germany), according to manufacturer's instructions. DNAse treatment (Qiagen, Hilden, Germany) was carried out to avoid possible contamination from genomic DNA. RNA purity and concentration were determined using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE) and RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). All samples had an RNA integrity number above 8.

**Library preparation and RNA-sequencing.** Eleven tagged libraries for RNA sequencing experiments were prepared from extracted mRNA using Agilent's SureSelect Strand Specific RNA Library Preparation Kit (Agilent Technologies), following the manufacture's instruction. Briefly, poly(A) mRNA was purified from 200 ng of total extracted RNA and fragmented using an RNA-seq Fragmentation Mix. First-strand and second-strand cDNA were synthetized and
end repaired. Adenylation of cDNA 3' ends and adaptor ligation were performed. Sixteen cycles of PCR amplification were used to amplify and index the adaptor-ligated cDNA library, and the PCR products were purified and size selected using the SPRIselect reagent kit (Beckman Coulter, Brea, CA, USA). A Qubit RNA Assay kit in a Qubit 2.0 Fluorometer (Life Technologies, CA, USA) was used to measure library concentrations. Individual libraries were monitored for insert size using the Agilent High Sensitivity DNA assay on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA). Multiplexed paired-end sequencing 2 × 100 bp was carried out on an Illumina Hi-Seq 2500 (UC Davis Genome Center, Davis, CA, USA). To this end, a total of 3 lanes were run. In order to avoid potential biases related to divergent sequencing qualities, four barcoded libraries were loaded in each lane by randomizing samples belonging to the three diets. Raw Illumina sequencing data have been deposited in GenBank (SRA) (PRJNA303144, BioProject SRP066556, study accession).

**Reads processing and quality control.** Initial quality control was carried out with the FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) version 0.11.4 [62]. Raw reads were then trimmed for low quality bases using CLC Genomics Workbench v 8.5 as follows: 1) all reads were trimmed of 10 nucleotides at the 5' end; 2) Illumina adapters were removed; 3) reads with >5% nucleotide with PHRED scores <20 were filtered out.

**Reads mapping.** The trimmed reads were mapped to the *Bos taurus* reference genome (UMD3.1, release 81.3) using CLC Genomics Workbench v 8.5 (CLC Bio, Aarhus, Denmark) with default parameters, except for a length fraction of 0.8 and a similarity fraction of 0.8. The alignments were saved as BAM files. In order to validate the steadiness of counts predicted by the CLC Genomics Workbench, gene count estimation was also performed using the HTSeq-count tool (part of the HTSeq framework, version 0.7.1, in 'union' mode), an open source toolkit that allows the input of raw counts from aligned reads to be annotated with gene names based on genomic features [63]. To compare the two methods, Pearson and non-parametric Spearman rank correlations were applied using the rcorr function in the Hmisc package in the R environment (v 3.1).

**Identification of alternative splicing events.** In order to improve the alternative splicing analysis, instead of using the previously obtained alignments, trimmed reads were mapped against the *Bos taurus* Ensembl genome UMD3.1 release 81.3 in a splicing-aware manner using STAR aligner v2.4.1d with the following arguments "-outSAMtype BAM SortedBy-Coordinate-outSAMmapqUnique 60-outFilterMultimapNmax 2-outFilterMismatchNmax 10-alignEndsType EndToEnd" [64]. Alternative 3' (A3SS) and 5' splice sites (A5SS), skipped exons (SE), mutually exclusive exons (MXE), and retained introns (RI) from the alignments were quantified using rMATS v3.2.0.beta [65] with Ensembl cow UMD3.1 annotation release 81.3. Fragment insert size and standard deviations, required to run rMATS, were calculated for each sample using the Picard CollectionInsertSizeMetric tool. Alternative splicing (AS) events were detected using both junction counts and reads on target. To identify different putative AS events due to the three different diets, rMATS evaluated whether differences in the exon or intron inclusion levels (IncLevelDifference) of an AS event between two conditions exceeds a false discovery rate (FDR) threshold of 5%.

**Identification of the differentially expressed genes.** The read counts obtained from the CLC Genomics Workbench were used to estimate exon expression and identify DEGs, which was done using the Bioconductor package edgeR version 3.10.0 [66] in the R environment (version 3.1). The "calcNormFactors" normalization function of the edgeR package was used to find a set of scaling factors for the library sizes that minimized the log-fold changes between samples. The scale factors were computed using a trimmed mean of M-values (TMM) between samples [66]. Principal component analysis (PCA) and Multidimensional Scaling analysis (MDS) were performed with FactoMineR and Ade4 packages (version 1.23) to assess whether
any particular array made a large contribution to the variability in the gene expression data, that is, whether it retained most of the information. The differences that emerged from the PCA analysis were assessed using the Monte Carlo Permutation Procedure (999 replicates; "randtest" function) of the Ade4 package in R [67]. DEGs were detected by applying a generalized linear model (GLM) likelihood ratio test, which fits negative binomial (NB) with the Cox-Reid dispersion estimates [67]. Common, trended and tagwise dispersions were calculated using the "estimateDisp" functions in the edgeR package. The DEGs were then determined using the GLM likelihood ratio test then each of the treatment groups was compared using the contrast argument of the "glmLRT" function in the edgeR package. EdgeR generated a list of genes with fold changes (FC), P-values, and the associated Benjamin-Hochberg FDR values. A further filtering step was introduced to delete transcripts with low expression levels (log2 < 4 in 90% of samples). The false discovery rate was then recalculated on the basis of the final transcript number, by implementing in R studio the FDR calculation method described by Storey and Tibshirani [68]. Only transcripts with FDR < 0.05 were retained for subsequent analyses. Visualization of the DEGs found to be commonly regulated in the three diet groups was performed using Venny, an interactive tool available at http://bioinfogp.cnb.csic.es/tools/venny/.

**Functional analysis of differentially expressed genes.** To gain a better understanding of the functional implications of these DEGs among the different dietary treatments, GO [69] and KEGG [70] enrichment analyses were performed using Cytoscape V2.7 (http://cytoscape.org/) with the ClueGo V1.3 plug-in 3 [71]. A single cluster analysis was used for each comparison. ClueGO determines the distribution of the gene list for the various GO terms and pathways. The P-value was calculated using right-sided hypergeometric tests and the Benjamini-Hochberg correction for multiple testing (FDR < 0.05). Together with this stringent FDR threshold, a high kappa value (0.4) enabled us to precisely select enriched GO terms in highly connected genes. The size of the nodes reflects the degree of enrichment of the terms. The network was automatically laid out using the organic layout algorithm in Cytoscape. Functional groups were created by iterative merging of the initially defined groups, according to the predefined kappa threshold. Only functional groups represented by their most significant term were visualized in the network [71]. A complementary approach, which involved feeding the list of selected genes into an IPA (ver. 5.5, Ingenuity Systems, Redwood City, CA), was also taken to identify relevant categories of molecular functions, cellular components and biological processes. Using this approach, we identified statistically overrepresented functional GO annotations, and determined their up- or down-expression and group-specific transcriptional networks. All listed or reconstructed cellular pathways were derived from the expert annotated database provided by the Ingenuity Knowledge Base. The IPA annotations follow the GO annotation principle, but are based on a proprietary knowledge database of over 10⁶ protein-protein interactions. The IPA output included biological functions and signaling pathways with statistical assessment of the significance of their representation based on Fisher's Exact test. IPA computed the networks and ranked them according to statistical likelihood [72]. Only the canonical pathways that presented a log2 (P-value) exceeding 1.30 (FDR < 0.05) were considered.

## Correlation analyses between differentially expressed genes and phenotypic traits

To find possible associations between the set of the previously identified DEG and phenotypic traits, the log2 of normalized gene counts for all DEGs expressed among the three dietary treatments were correlated with: i) animal growth performance, ii) carcass traits, iii) liver weight and chemical composition, iv) liver mineral composition, and v) liver FA composition.
Pearson product-moment correlation coefficients and *P*-values were calculated using the Hmisc package (3.17.2) in the R environment. The level of significance was set at *P* < 0.05. Heatmaps and dendograms were generated using the significant correlations for at least one trait by the gplots package in implemented in the R environment (v 3.1).

**Supporting Information**

**S1 Fig.** Boxplot of the natural log transformed estimated gene counts per animal after trimmed mean of M-values (TMM) normalization. Boxplots were created with the ggplot2 package in the R statistical environment.
(TIFF)

**S2 Fig.** Interclass principal component analysis (PCA) of normalized gene count data. The plot was created using principal component analysis function implemented in the ade4 R package with samples as variables and differential expression comparison groups as class levels.
(TIFF)

**S1 Table.** Mapping summary statistics of RNA sequencing data. Comparative results obtained by CLC Genomics Workbench and HTSeq-count tools.
(XLS)

**S2 Table.** Differential alternative splicing events for all comparisons using rMATS. Significant differences (FDR < 0.05) in the inclusion of splicing events across the experimental diets. Conv: conventional diet; LowI: low-impact diet.
(XLSX)

**S3 Table.** Differentially expressed genes (FDR < 0.05) for comparison among experimental diets. LowI: low-impact; Conv: conventional diet; Nutr: nutraceutical diet.
(XLSX)

**S4 Table.** Canonical pathways enriched in the experimental diets using IPA. The ratio values (number of molecules in a given pathway that meets cut criteria, divided by total number of molecules that make up that pathway) are also presented.
(XLSX)

**S5 Table.** Upstream regulators analysis results for comparison among experimental diets using IPA. IPA prediction of upstream regulators that are activated or inhibited to explain the upregulated and downregulated genes observed in the experimental diets.
(XLSX)

**S6 Table.** Correlation matrix of Pearson product-moment correlations between liver mineral profiles and differentially expressed genes. Pearson product-moment correlation coefficients and *P* values are reported for all correlations.
(TXT)

**Acknowledgments**

This research was supported by the University of Padova (Ricerca Scientifica fondi quota EX 60% -649 60A08-7943/14).

**Author Contributions**

**Conceptualization:** AC SS GB.
**Data curation:** SP MP MB.

**Formal analysis:** SP MP MB NM AC.

**Funding acquisition:** GB.

**Investigation:** SP AC MP NM.

**Methodology:** SP NM.

**Project administration:** GB AC SS.

**Resources:** GB SS LB.

**Software:** SP MP MB NM.

**Supervision:** GB AC SS.

**Visualization:** SP GB AC SS NM MP MB LB.

**Writing – original draft:** SP NM MP.

**Writing – review & editing:** SP GB AC SS NM MP MB LB.

## References

1. Horizon 2020—European Commission. https://ec.europa.eu/programmes/horizon2020/.

2. Schiavon S, Tagliapietra F, Dalla Montà G, Cecchinato A, Bittante G. Low protein diets and rumen-protected conjugated linoleic acid increase nitrogen efficiency and reduce the environmental impact of double-muscled young Piemontese bulls. Anim Feed Sci Technol. 2012; 174:96–107.

3. Schiavon S, Tagliapietra F, Cesaro G, Gallo L, Cecchinato A, Bittante G. Low crude protein diets and phase feeding for double-muscled crossbred young bulls and heifers. Livest Sci. 2013; 157:462–470.

4. López-Alonso M: Trace minerals and livestock: not too much not too little. ISRN Vet Sci. 2012; 2012:704825. doi: 10.5402/2012/704825 PMID: 23762589

5. Kris-Etherton PM. Fish Consumption, Fish Oil, Omega-3 Fatty Acids, and Cardiovascular Disease. Circulation. 2002; 106:2747–2757. PMID: 12438303

6. Mattos R: Effects of dietary fatty acids on reproduction in ruminants. Rev Reprod 2000, 5:38–45. PMID: 10711734

7. Mir Z, Rushfeldt ML, Mir PS, Paterson LJ, Weselake RJ. Effect of dietary supplementation with either conjugated linoleic acid (CLA) or linoleic acid rich oil on the CLA content of lamb tissues. Small Rumin Res. 2000; 36:25–31.

8. Nuernberg K, Dannenberger D, Nuernberg G, Ender K, Voigt J, Scollan ND, et al. Effect of a grass-based and a concentrate feeding system on meat quality characteristics and fatty acid composition of longissimus muscle in different cattle breeds. Livest Prod Sci. 2005; 94: 137–147.

9. Mach N, Zom RLG, Widjaja HCA, van Wikselaar PG, Weurding RE, Goselink RMA, et al. Dietary effects of linseed on fatty acid composition of milk and on liver, adipose and mammary gland metabolism of periparturient dairy cows. J Anim Physiol Anim Nutr (Berl). 2013; 97:89–104.

10. Corl B, Barbano DM, Bauman DE, Ip C. cis-9, trans-11 CLA derived endogenously from trans-11 18:1 reduces cancer risk in rats. J Nutr. 2003; 133:2893–2900. PMID: 12949384

11. Park Y, Storkson JM, Albright KJ, Liu W, Pariza MW. Evidence that the trans-10,cis-12 isomer of conjugated linoleic acid induces body composition changes in mice. Lipids. 1999; 34:235–41. PMID: 10230716

12. Azain MJ, Hausman DB, Sisk MB, Flatt WP, Jewell DE. Dietary conjugated linoleic acid reduces rat adipose tissue cell size rather than cell number. J Nutr. 2000; 130:1548–1554. PMID: 10827208

13. Baumgard LH, Matitashvili E, Corl BA, Dwyer DA, Bauman DE. trans-10, cis-12 conjugated linoleic acid decreases lipogenic rates and expression of genes involved in milk lipid synthesis in dairy cows. J Dairy Sci. 2002; 85:2155–63. doi: 10.3168/jds.S0022-0302(02)74294-X PMID: 12362447

14. Coyne GS, Kenny DA, Childs S, Sreenan JM, Waters SM. Dietary n-3 polyunsaturated fatty acids alter the expression of genes involved in prostaglandin biosynthesis in the bovine uterus. Theriogenology. 2008; 70:772–782. doi: 10.1016/j.theriogenology.2008.05.048 PMID: 18582926
15. Ahnadi CE, Beswick N, Delbecchi L, Kennelly JJ, Lacasse P. Addition of fish oil to diets for dairy cows. II. Effects on milk fat and gene expression of mammary lipogenic enzymes. J Dairy Res. 2002; 69: 521-531. PMID: 12463690

16. Baldwin RL, Li RW, Li C, Thomson JM, Bequette BJ. Characterization of the longissimus lumborum transcriptome response to adding propionate to the diet of growing Angus beef steers. Physiol Genomics. 2012; 44:543-50. doi: 10.1152/physiolgenomics.00144.2011 PMID: 22454452

17. Baldwin RL VI, Wu S, Li W, Li C, Bequette BJ, Li RW. Quantification of transcriptome responses of the rumen epithelium to butyrate infusion using RNA-seq technology. Gene Regul Syst Bio. 2012; 6: 67-80. doi: 10.4137/GRSB.S9687 PMID: 22654504

18. Matthews JC, Zhang Z, Patterson JD, Bridges PJ, Stromberg AJ, Boling JA. Hepatic transcriptome profiles differ among maturing beef heifers supplemented with inorganic, organic, or mixed (50% inorganic:50% organic) forms of dietary selenium. Biol Trace Elem Res. 2014; 160: 321-339. doi: 10.1007/s12011-014-0050-4 PMID: 24996959

19. da Costa AS, Bessa RJ, Pires VM, Rolo EA, Pinto RM, Andrade Fontes CM, et al. Is hepatic lipid metabolism of beef cattle influenced by breed and dietary silage level? BMC Vet Res. 2014; 10:65. doi: 10.1186/1746-6148-10-65 PMID: 24621212

20. Gruffat D, Gobert M, Durand D, Bauchart D. Distinct metabolism of linoleic and linolenic acids in liver and adipose tissues of finishing Normande cull cows. Animal. 2011; 5:1090-1098. doi: 10.1017/S1751731111000073 PMID: 22440104

21. Schiavon S, Tagliapietra F, Maso MD, Bailoni L, Bittante G. Effects of low-protein diets and rumen-protected conjugated linoleic acid on production and carcass traits of growing double-muscled Piemontese bull. J Anim Sci. 2010; 88:3372-3383 doi: 10.2527/jas.2009-2558 PMID: 20601521

22. Gleeson M, Bishop NC. Modification of immune responses to exercise by carbohydrate, glutamine and anti-oxidant supplements. Immunol Cell Biol. 2000; 78:554-561. doi: 10.1111/j.1440-1711.2000.t01-6-x PMID: 11050539

23. Berg JM, Tymoczko JL, Stryer L. Food Intake and Starvation Induce Metabolic Changes. In: Biochemistry, 5th edition. 2002.

24. Roden M, Perseghin G, Petersen KF, Hwang JH, Cline GW, Gerow K, et al. The roles of insulin and glucagon in the regulation of hepatic glycogen synthesis and turnover in humans. J Clin Invest. 1996; 97:642-8. doi: 10.1172/JCI118460 PMID: 8609218

25. Schiavon S, De Marchi M, Tagliapietra F, Bailoni L, Cecchinato A, Bittante G. Effect of high or low protein ration combined or not with rumen protected conjugated linoleic acid (CLA) on meat CLA content and quality traits of double-muscled Piemontese bulls. Meat Sci. 2011; 89:133-142. doi: 10.1016/j.meatsci.2011.03.025 PMID: 21561723

26. Madron MS, Peterson DG, Dwyer DA, Corl BA, Baumgard LH, Beermann DH, et al. Effect of extruded full-fat soybeans on conjugated linoleic acid content of intramuscular, intermuscular, and subcutaneous fat in beef steers. J Anim Sci. 2002; 80:1135-1143. PMID: 12002322

27. Pellattiero E, Cecchinato A, Tagliapietra F, Schiavon S, Bittante G. Determination by GCxGC of fatty acid and conjugated linoleic acid (CLA) isomer profiles in six selected tissues of lambs fed on pasture or on indoor diets with and without rumen-protected CLA. J Agric Food Chem. 2015; 63:963-974. doi: 10.1021/jf504956x PMID: 25525905

28. Bernard L, Leroux C, Chilliard Y. Expression and nutritional regulation of lipogenic genes in the ruminant lactating mammary gland. Adv Exp Med Biol. 2008; 606: 67-108. doi: 10.1007/978-0-387-74087-4_2 PMID: 18183925

29. López Alonso M, Benedito JL, Miranda M, Castillo C, Hernández J, Shore RF. Toxic and trace elements in liver, kidney and meat from cattle slaughtered in Galicia (NW Spain). Food Addit Contam. 2000; 17:447-57. doi: 10.1080/02652030050034028 PMID: 10932787

30. Alonso ML, Montana FP, Miranda M, Castillo C, Hernandez J, Benedito JL. Interactions between toxic (As, Cd, Hg and Pb) and nutritional essential (Ca, Co, Cr, Cu, Fe, Mn, Mo, Ni, Se, Zn) elements in the tissues of cattle from NW Spain. BioMetals. 2004; 17: 389-397. PMID: 15259359

31. Kotula AW, Lusby WR. Mineral composition of muscles of 1- to 6-year-old steers. J Anim Sci. 1982; 54: 544-8. PMID: 7085514

32. Kreipe L, Vernay MCMB, Oppliger A, Wellnitz O, Bruckmaier RM, van Dorland HA. Induced hypoglycemia for 48 hours indicates differential glucose and insulin effects on liver metabolism in dairy cows. J Dairy Sci. 2011; 94: 5435-5448. doi: 10.3168/jds.2011-4208 PMID: 22032366

33. Davidson AJ, Castañón-Cervantes O, Stephan FK. Daily oscillations in liver function: diurnal vs circadian rhythmicity. Liver Int. 2004; 24:179-186. doi: 10.1111/j.1478-3231.2004.0917.x. PMID: 15189266

34. Fernandez ML, West KL. Mechanisms by which dietary fatty acids modulate plasma lipids. J Nutr. 2005; 135:2075-2078. PMID: 16140878
35. Calder PC. n-3 Polyunsaturated fatty acids, inflammation, and inflammatory diseases. Am J Clin Nutr. 2006; 83:S1505-1519.

36. Cannell IG, Kong YW, Bushell M. How do microRNAs regulate gene expression? Biochem Soc Trans. 2008; 36(Pt 6):1224-31. doi: 10.1042/BST0361224 PMID: 19021530

37. Yang M, Liu W, Pellicane C, Sahyoun C, Joseph BK, Gallo-Ebert C, Donigan M, et al. Identification of miR-185 as a regulator of de novo cholesterol biosynthesis and low density lipoprotein uptake. J Lipid Res. 2014; 55:226-238. doi: 10.1194/jlr.M041335 PMID: 24296663

38. Dong X-Y, Tang S-Q, Chen J-D. Dual functions of Insig proteins in cholesterol homeostasis. Lipids Health Dis. 2012; 11:173. doi: 10.1186/1476-511X-11-173 PMID: 23249523

39. Loewen CJ., Levine TP. Cholesterol Homeostasis: Not until the SCAP Lady INSIGs. Curr Biol. 2002; 12:R779-R781. PMID: 12445404

40. Jacobs AAA, Dijkstra J, Liesman JS, Vandehaar MJ, Lock AL, van Vuuren AM, et al. Effects of short- and long-chain fatty acids on the expression of stearoyl-CoA desaturase and other lipogenic genes in bovine mammary epithelial cells. Animal. 2013; 7:1508-16. doi: 10.1017/S175173111300061X PMID: 23597233

41. Mach N, Jacobs AA, Kruijst L, van Baal J, Smits MA. Alteration of gene expression in mammary gland tissue of dairy cows in response to dietary unsaturated fatty acids. Animal. 2011; 5:1217-30. doi: 10.1017/S1751731111000103 PMID: 22440174

42. Hiller B, Angulo J, Olivera M, Nuernberg G, Nuernberg K. How selected tissues of lactating holstein cows respond to dietary polyunsaturated fatty acid supplementation. Lipids. 2013; 48:357-67. doi: 10.1007/s11745-012-3737-3 PMID: 23129256

43. Matsunami T, Sato Y, Ariga S, Sato T, Shimomura T, Kashimura H, et al. Regulation of synthesis and oxidation of fatty acids by adiponectin receptors (AdipoR1/R2) and insulin receptor substrate isoforms (IRS-1/-2) of the liver in a nonalcoholic steatohepatitis animal model. Metabolism. 2011; 60:805-14. doi: 10.1016/j.metabol.2010.07.032 PMID: 20846698

44. Jeoung NH, Wu P, Joshi MA, Jaskiewicz J, Bock CB, Depaoli-Roach AA, et al. Role of pyruvate dehydrogenase kinase isoenzyme 4 (PDHK4) in glucose homeostasis during starvation. Biochem J. 2006; 397:417-25. doi: 10.1042/BJ20060125 PMID: 16606348

45. Jump DB.N-3 polyunsaturated fatty acid regulation of hepatic gene transcription. Curr Opin Lipidol. 2008; 19:242-7. doi: 10.1097/MOL.0b013e3282ffaf6a PMID: 18460914

46. Schmitz G, Ecker J. The opposing effects of n-3 and n-6 fatty acids. Prog Lipid Res. 2008; 47:147-55. doi: 10.1016/j.plipres.2007.12.004 PMID: 18198131

47. Frost RJA, Olson EN. Control of glucose homeostasis and insulin sensitivity by the Let-7 family of micro-RNAs. Proc Natl Acad Sci. 2011; 108:21075-21080. doi: 10.1073/pnas.1118922109 PMID: 22160727

48. Sun C, Huang F, Liu X, Xiao X, Yang M, Hu G, et al. miR-21 regulates triglyceride and cholesterol metabolism in non-alcoholic fatty liver disease by targeting HMGCR. Int J Mol Med. 2015; 35:847-53. doi: 10.3892/ijmm.2015.2076 PMID: 25605429

49. Ravi S, Schilder RJ, Kimball SR. Role of Precursor mRNA Splicing in Nutrient-Induced Alterations in Gene Expression and Metabolism. J Nutr. 2015; 145:841-846. doi: 10.3945/jn.114.203216 PMID: 25761502

50. Okamura Y, Watari M, Jerud ES, Young DW, Ishizaka ST, Rose J, et al. The Extra Domain A of Fibronectin Activates Toll-like Receptor 4. J Biol Chem. 2001; 276:10229-10233. doi: 10.1074/jbc.M100099200 PMID: 11150311

51. Bagashev A, Fitzgerald MC, LaRosa DF, Rose PP, Cherry S, Johnson AC, et al. Leucine-Rich Repeat (in Flightless I) Interacting Protein-1 Regulates a Rapid Type I Interferon Response. J Interferon Cyto-kine Res. 2010; 30: 843-852. doi: 10.1089/jir.2010.0017 PMID: 20586614

52. Markiewski MM, Lambris JD. The role of complement in inflammatory diseases from behind the scenes into the spotlight. Am J Pathol. 2007; 171: 715-27. doi: 10.2353/ajpath.2007.070166 PMID: 17640961

53. White ES, Muro AF. Fibronectin splice variants. Understanding their multiple roles in health and disease using engineered mouse models. IUBMB Life. 2011; 63:538-546. doi: 10.1002/iub.493 PMID: 21698758

54. Cildir G, Akıncılar SC, Tergaonkar V. Chronic adipose tissue inflammation: all immune cells on the stage. Trends Mol Med. 2013; 19:487-500. doi: 10.1016/j.molmed.2013.05.001 PMID: 23746697

55. Yang X, Downes M, Yu RT, Bookout AL, He W, Straume M, et al. Nuclear receptor expression links the circadian clock to metabolism. Cell. 2006; 126:801-10. doi: 10.1016/j.cell.2006.06.050 PMID: 16923398

56. Boyer JL, Graf J, Meier PJ. Hepatic transport systems regulating pHi, cell volume, and bile secretion. Annu Rev Physiol. 1992; 54:415-38. doi: 10.1146/annurev.ph.54.030192.002215 PMID: 1562181
57. Palmiter RD. Protection against zinc toxicity by metallothionein and zinc transporter 1. Proc Natl Acad Sci U S A. 2004; 101:4918-4923. doi: 10.1073/pnas.0401022101 PMID: 15041749

58. Dacasto M, Eeckhoutte C, Capolongo F, Dupuy J, Carletti M, Calléja C, et al. Effect of breed and gender on bovine liver cytochrome P450 3A (CYP3A) expression and inter-species comparison with other domestic ruminants. Vet Res. 2005; 36:179-90. doi: 10.1051/vetres:2004066 PMID: 15720971

59. Nakamura TJ, Sellix MT, Kudo T, Nakao N, Yoshimura T, Ebihara S, et al. Influence of the estrous cycle on clock gene expression in reproductive tissues: effects of fluctuating ovarian steroid hormone levels. Steroids. 2010; 75:203-12.2. doi: 10.1016/j.steroids.2010.01.007 PMID: 20096720

60. Jenkins TC. Technical note: common analytical errors yielding inaccurate results during analysis of fatty acids in feed and digesta samples. J Dairy Sci. 2010; 93:1170-4. doi: 10.3168/jds.2009-2509 PMID: 20172237

61. Schiavon S, Pellattiero E, Cecchinato A, Tagliapietra F, Dannenberger D, Nuernberg K, et al. Influence of different sample preparation procedures on the gas chromatographic fatty acid profile of subcutaneous fat, liver and muscle of beef cattle. J Dairy Sci. 2016; Accepted manuscript.

62. Andrews, S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/.

63. Anders S, McCarthy DJ, Chen Y, Okoniewski M, Smyth GK, Huber W, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013; 8:1765-1786. doi: 10.1038/nprot.2013.099 PMID: 23975260

64. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013; 29:15-21. doi: 10.1093/bioinformatics/bts635 PMID: 23104886

65. Shen S, Park JW, Lu Z, Lin L, Henry MD, Wu YN, et al. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci U S A. 2014; 111:E5593-601. doi: 10.1073/pnas.1419161111 PMID: 25480548

66. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26:139-40. doi: 10.1093/bioinformatics/btp616 PMID: 19910308

67. Dray S, Chessel D, Thioulouse J. Co-inertia analysis and the linking of ecological data tables. Ecology. 2003; 84: 3078-3089.

68. Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A. 2003; 100:9440-5. doi: 10.1073/pnas.1530509100 PMID: 12883005

69. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000; 25:25-9. doi: 10.1038/75556 PMID: 10802651

70. Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 1999; 27:29-34. PMID: 9847135

71. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009; 25:1091-3. doi: 10.1093/bioinformatics/btp101 PMID: 19237447

72. Calvano SE, Xiao W, Richards DR, Felciano RM, Baker H V., Cho RJ, et al. A network-based analysis of systemic inflammation in humans. Nature. 2005; 437:1032-1037. doi: 10.1038/nature03985 PMID: 16136080