Effect of a diet enriched with omega-6 and omega-3 fatty acids on the pig liver transcriptome

Abstract

The optimal ratio of omega-6 to omega-3 polyunsaturated fatty acids (PUFAs) is important for keeping the homeostasis of biological processes and metabolism, yet the underlying biological mechanism is poorly understood. The objective of this study was to identify changes in the pig liver transcriptome induced by a diet enriched with omega-6 and omega-3 fatty acids and to characterize the biological mechanisms related to PUFA metabolism.

Polish Landrace pigs (n = 12) were fed diet enriched with linoleic acid (LA, omega-6) and α-linolenic acid (ALA, omega-3) or standard diet as a control. The fatty acid profiling was assayed in order to verify how feeding influenced the fatty acid content in the liver, and subsequently next-generation sequencing (NGS) was used to identify differentially expressed genes (DEG) between transcriptomes between dietary groups. The biological mechanisms and pathway interaction networks were identified using DAVID and Cytoscape tools. Fatty acid profile analysis indicated a higher contribution of PUFAs in the liver for LA- and ALA-enriched diet group, particularly for the omega-3 fatty acid family, but not omega-6. Next-generation sequencing identified 3565 DEG, 1484 of which were induced and 2081 were suppressed by PUFA supplementation. A low ratio of omega-6/omega-3 fatty acids resulted in the modulation of fatty acid metabolism pathways and over-representation of genes involved in energy metabolism, signal transduction, and immune response pathways.

In conclusion, a diet enriched with omega-6 and omega-3 fatty acids altered the transcriptomic profile of the pig liver and would influence animal health status.

Introduction

A balanced ratio of omega-6 to omega-3 polyunsaturated fatty acids (PUFAs) plays an essential role in multiple biological processes and maintains the metabolic homeostasis. It is widely accepted that long-chain PUFAs are considered as a bioactive molecules, which serve beneficial outcomes to health. However, the mechanism responsible for a health-promoting effect as a result of fatty acid consumption is still not fully clarified.

Polyunsaturated fatty acids, particularly omega-6 or omega-3 family, counteract cardiovascular (Stanley et al. 2013) and metabolic disorders (Warensjo et al. 2006). Efficacy of omega-3-rich diet for reducing blood cholesterol levels and decreasing risk of cardiovascular disease, arthritis, neurodegenerative diseases (Lands 2012), and risk of progression to psychotic disorders (Amminger et al. 2015) is unambiguous. Omega-3 PUFAs derivatives can act as signaling molecules, modulating the anti-inflammatory response, and controlling cellular processes involved in programmed cell death (apoptosis), lymphocyte proliferation, suppression of inflammatory cytokines, and phagocytosis (de Pablo and Alvarez de Cienfuegos 2000; Kruger et al. 2010; Wymann and Schneiter 2008). Since omega-6 PUFAs are involved in the regulation of eicosanoid synthesis, they control the activity of the immune system. They are also associated with a reduced risk of cardiovascular diseases. Many recent studies have investigated fish oils and plant oils (for instance linseed and rapeseed oils) rich in omega-3 fatty acids, indicating the protective actions against the development of cardiovascular and metabolic diseases such as obesity, type II diabetes, and metabolic syndrome (Hanke et al. 2013; Lionetti et al. 2014; Morine et al. 2011). Especially, a low omega-6/omega-3 ratio decreases the risk of metabolic diseases (Jump 2011).

In regard to plant oils, diet can provide short-chain linoleic acid (LA) and α-linolenic acid (ALA), which are the main precursors of omega-6 and omega-3 fatty acids, respectively, and are considered as essential fatty acids for humans and other mammals. They can be stored in tissues or converted into long-chain fatty acids such as arachidonic acid (AA), eicosapentaenoic acid (EPA), and docosahexaenoic acid (DHA), which are subsequently involved in the synthesis of hormones such as prostaglandins, prostacyclins, thromboxanes, or leukotrienes (Aro 2003; Harris et al. 2006). There is now significant evidence that one of the mechanisms by which omega-6 and omega-3 fatty acids exert their protective role is interaction with nuclear receptors. They are the natural ligands that influence gene expression and regulate the activity of numerous transcription factors, including sterol regulatory element-binding protein (SREBP), peroxisome proliferator-activated receptor (PPAR), liver X receptor (LXR), farnesoid X receptor (FXR), and retinoid X receptor (RXR) (Getek et al. 2013; Gormaz et al. 2010; Vallim and Salter 2010).

In our studies, we have chosen the pig as an animal model for humans. Such a model has additional advantage because of the relatively similar anatomy, body morphology, physiological and metabolic processes to humans as well as high sequence homology, and a similar chromosomal structure to the human genome, compared to other domestic animal species (Lunney 2007). Therefore, processes revealed in the pig transcriptome are more likely to be similar to those in humans. For these reasons, the pig is considered to be a good animal model for nutrigenomic and metabolic research.

Since the liver is an organ which plays a major role in the regulation of the availability of lipids for other tissues by modulating the uptake and synthesis of lipids and the secretion of lipoproteins, it is adequate to study the fatty acid metabolism dysregulations (Nguyen et al. 2008). Besides adipose tissue, the liver is one of the main locations of synthesis for both cholesterol and triglyceride fatty acids (Bergen and Mersmann 2005). The synthesis of triglycerides depends largely on the availability of fatty acids, which are derived from de novo synthesis, received from diet, or released from tissues (Vallim and Salter 2010). Fatty acid oxidation occurs in hepatic mitochondria and peroxisomes. Dietary fatty acids affect de novo lipogenesis in the liver by controlling the expression of the key enzymes of fatty acid synthesis, particularly acetyl-CoA carboxylase (ACC) and fatty acid synthase (FAS) by PUFA (Guillevic et al. 2009) or stearoyl-CoA desaturase (SCD) by monounsaturated fatty acids (MUFA) (Ntambi 1999).

Nowadays, high-throughput technology, such as next-generation sequencing (NGS), becomes increasingly popular, providing identification of changes in the transcriptome in nutrigenomics investigations (Chen et al. 2011; Li et al. 2011), further validated by standard quantitative real-time PCR (qRT-PCR) (Nygard et al. 2007; Shi and Chiang 2005). The objective of this study was the use of the RNA-Seq approach to identify changes in the pig liver transcriptome induced by a diet enriched with LA (omega-6 family) and ALA (omega-3) fatty acids and to indicate the general biological mechanisms underlying the effects related to PUFA metabolism.

Materials and methods

Animals and treatments

The experiment was carried out on Polish Landrace pigs (n = 12) housed at a commercial farm (Wronie, Poland) in two separate pens under standard conditions. The first stage of feeding was composed of standard diet appropriate for age group and included all of the pigs, which were fed until reaching 60 kg in body weight (BW) (see scheme in Fig. 1). In the second fattening period, the pigs were divided into two feeding groups (six pigs per group) and fed diets differing in omega-6 and omega-3 fatty acid content based on linseed and rapeseed oil supplementation. In the control group (C, n = 6), feed contained 268 mg LA/100 g (linoleic acid—omega-6) and 25 mg ALA/100 g (α-linolenic acid—omega-3 family). The experimental diet (E, n = 6) was enriched with LA and ALA and contained 660 mg and 64 mg/100 g diet, respectively. The feeding experiment was performed on an equal number of pigs from each treatment group, and all experiment included gilts. Diets were provided ad libitum with a constant access to water. The fodder mixtures were isoenergetic (13 MJ EM per kg of mixture) and isoprotein balanced (15.5 % crude protein), while differing in fat content ranging from 1.78 % (C group) to 4.72 % (E group). Diets were composed of the same components, except linseed and rapeseed oils, which were provided as a polyunsaturated fatty acid source (Table 1). The feed was administered for approximately 4 months until slaughter at the average BW of 110 kg. All animals were characterized by aligned age and body weight. Pigs were sacrificed by electrical stunning and exsanguination, according to industry standards. After slaughter, liver samples were excised and immediately frozen in liquid nitrogen and further stored at −80 °C. Slaughter procedure was carried out with the required permits and according to Minister of Agriculture and Rural Development dated April 2, 2004.

Fig. 1

Experimental design. The first stage of feeding was composed of a standard diet until the pigs reached 60 kg in body weight. In the second stage, the pigs were divided into two groups (six pigs per group): the experimental (E) group—fed a diet enriched with LA and ALA, which contained 660 mg and 64 mg/100 g diet, respectively; and the control (C) group—fed a standard diet which contained 268 mg LA/100 g and 25 mg ALA/100 g

Table 1 Raw material composition of mixtures (%)

Fatty acid profiling

The fatty acid composition in the liver was assayed in order to verify how feeding influenced the fatty acid content in tissues. Subsequently, the ratio of omega-6/omega-3 fatty acids in the liver was calculated for both dietary groups and used as the basis for further investigation. For p-value determination, unpaired Student’s test was used.

The fatty acid profile analysis was performed using the chloroform-methanol (2:1 v/v) method (Folch et al. 1957) on 1 g of liver sample. Fatty acid methyl esters were determined by gas chromatography (GC-FID) on the Agilent 7890 Series GC system equipment with a Hewlett-Packard-88 capillary column (60 m × 0.25 mm × 0.20 μm) (Agilent, J&W GC Columns, USA) with helium as a carrier gas. The details of the analysis are described in Polawska et al. (2013). Fatty acids were identified by comparison with fatty acid standards (Supelco 37 Component FAME mix, 47885-U—10 mg/ml in methylene chloride, analytical standard, Sigma-Aldrich Co.) and expressed in g/100 g fatty acid methyl esters (FAME).

Gene expression analysis

Next-generation sequencing

The total RNA for sequencing were purified from 20 mg of tissue sample of 12 pigs (6 animals in control group and 6 animals in PUFA treatment group) using the messenger RNA (mRNA) extraction kit RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s recommendations. The concentration and purity of total RNA were checked using the UV/Vis spectrophotometer (Actgene, New Jersey, USA) at 260 nm, and the integrity was evaluated using a microfluidic assay on Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). Only high quality RNA extracts (RNA integrity number (RIN) ≥8) were used for pooling within each treatment group using equal amounts of RNA per animal.

Next-generation sequencing data were obtained from pooled RNA samples within each group to ensure the most robust transcriptome. Construction of complementary DNA (cDNA) libraries for RNA sequencing was prepared using the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA). RNA-Seq analysis was performed for identification of transcriptional changes using MiSeq instrument (Illumina) according to the manufacturer’s recommendations at the laboratory of Genomed (http://www.genomed.pl/) using paired-end libraries. Two replicates of each pool independent for library synthesis and sequencing were done. The number of reads equivalent to mapped reads (reads per kilobase per million (RPKM)) was used to normalize the expression of each gene. Benjamini and Hochberg test (an error rate 0.05) and false discovery rate correction test (FDR < 0.05) were used for adjusting the p-value. In this study, a minimum twofold change (FC ≥ 2) difference in gene expression between treatment groups was used as the inclusion criteria for a gene. Reads were mapped to the Sus scrofa v. 10.2 genome as a reference genome. The quality of obtained data was checked in regard to the presence and abundance of contaminating sequences, average read length, and GC content. The NGS experiment conformed to the MIAME guidelines (Brazma et al. 2001) except for the microarray design. The RNA-Seq raw data were deposited in NCBI GEO (Gene Expression Omnibus) (Barrett et al. 2013) under accession number GSE72123.

Quantitative real-time PCR

Total cellular RNA was extracted according to the Chomczyński protocol (Chomczynski and Sacchi 1987; Chomczynski and Sacchi 2006). The samples of liver tissue (20 mg) were lysed with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and homogenized on MagNA Lyser (Roche Diagnostics, Mannheim, Germany). The concentration and purity of the RNA extracts were measured using UVS-99 micro-volume UV/Vis spectrophotometer (ACTgene, Piscataway, NJ, USA) at 260 nm. Integrity was checked using the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA) based on the RNA integrity number (RIN) parameter. The RIN values for all samples were higher than 7.5. Genomic DNA was removed from RNA extracts by incubation with RQ1 RNase-free DNase (Promega, Madison, WI, USA) for 15 min at 25 °C. The RNA extracts were stored in RNase-free water at −80 °C. The reverse transcription reaction was performed using the Transcriptor First Strand cDNA Synthesis Kit (Roche Applied Science, Penzberg, Germany) from 1.5 μg of total RNA in a final volume of 20 μl using the oligo-dT primers.

Quantitative real-time PCR was performed on 30 individual samples (biological replicates—15 samples for each dietary group) to validate the reliability of RNA-Seq gene expression data. For confirmation, we have chosen the differentially expressed genes EXTL1, APOA4, APOA5, ACSL1, ELOVL6, FASN, COL1A1, DGAT2, FADS1, FADS2, and MMP2. Most of the selected genes are involved in lipid metabolism. Gene expression was normalized by the use of topoisomerase (DNA) II beta (TOP2B) (Nygard et al. 2007; Pierzchala et al. 2011) as a most stable gene revealed by NormFinder algorithm, which was used to determine the optimal normalization gene among a set of candidates (GAPDH, TBP, TOP2B). The primer sets for real-time PCR were designed using Primer-BLAST tool (Ye et al. 2012) (www.ncbi.nlm.nih.gov/tools/primer-blast/) and Primer Premier 5.0 software (PREMIER Biosoft, Palo Alto, CA, USA). The sequences of the primers and their annealing temperatures are described in Table 2.

Table 2 Oligonucleotide primer sequences used for quantitative real-time PCR experiments

The reaction mixture for real-time PCR consisted of 4.5 μl of cDNA template (1/20 dilution), 5 pmol reaction primers, and 5 μl SYBR Green I in a final volume of 10 μl. All PCR reactions were performed in duplicate. For each gene, negative controls without the cDNA template were included. The reactions were performed in a 96-well optical plate on the Light Cycler 96 instrument (Roche Applied Science, Germany). Cycling conditions were as follows: pre-incubation at 95 °C for 5 min followed by 45 amplification cycles of 95 °C (10 s), 59 °C (10 s), and 72 °C (10 s). After each experiment, the qRT-PCR products were verified for specific amplification by melting curve analysis: 95 °C (5 s), 70 °C (1 min), heating to 95 °C, and cooling step to 37 °C for 30 s. Additionally, PCR product quality and specificity with regard to the applied primer pairs were checked using gel electrophoresis.

The gene expression levels were calculated using the LinRegPCR analysis program according to the baseline estimation method (Ramakers et al. 2003; Ruijter et al. 2013; Ruijter et al. 2009). In order to evaluate the relative gene expressions for independent samples, bilateral Student’s test was used and p-value <0.05 was considered significant. The quantitative real-time PCR was performed according to the MIQE guidelines (Bustin et al. 2013).

Bioinformatics analysis

For transcriptome assembling and dispersion evaluation, the CLC Genomics workbench v. 6.0. was used. Since pooled samples were used in the NGS experiment, the expression variance for each gene was estimated by using the variance across the four sequencing for that gene. This “blind” method of dispersion estimate is used in the absence of biological replicates, even though it can overly conservative statistical significance testing of differential gene expression (Anders and Huber 2010).

All data processing steps were performed using Lasergene Genomic Suite v. 5 software package (DNASTAR, Madison, WI, USA) for differential gene expression evaluation. The first step included selection of the genes that showed high variation for all identified gene transcripts. Only genes that demonstrated more than a twofold change in differential expression between dietary groups with a p-value ≤0.05 for the T test were selected for further analysis. Additionally, genes that demonstrated dietary group-specific expression (and remained below the detection level in the opposite dietary group) were investigated separately.

Biological mechanisms underlying DEGs were investigated using the DAVID v. 6.7 software (the Database for Annotation, Visualization, and Integrated Discovery) (http://david.abcc.ncifcrf.gov/) (Huang et al. 2009a; Huang et al. 2009b). The sets of genes were uploaded using ENTREZ gene IDs. To obtain maximal annotation results, the pig gene IDs were converted to the homologues human gene IDs. The p-values of the enrichment of the number of genes in biological mechanisms were evaluated by Benjamini correction, and values less than 0.05 were considered significant. The data were analyzed for the following situations: (1) all DEGs, (2) DEGs with expression levels below the detection level in the C group, and (3) DEGs with expression levels below the detection level in the E group. The analysis performed using the DAVID software allowed to derive the comprehensive information about the genes’ participation in biological pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) (Kanehisa and Goto 2000; Kanehisa et al. 2014). The p-values calculated in the DAVID software indicate the importance of the general biological mechanisms revealed in all animals but do not distinguish the inter-individual differences. For this purpose, we examined the alterations in the pig transcriptome in pooled liver samples per treatment group in order to clear out the bias caused by variation in a heterogeneous group of individuals derived from the general population. Apart from RNA-Seq, qRT-PCR was performed for individual samples, and indeed, high variation between individual animals was observed. The physiological effects of dietary omega-6 and omega-3 fatty acid intake were measured in the individual animals which are the consequence of the interactions between the animals’ genotype and the omega-6/omega-3 ratio. Therefore, it was expected that animal-specific mechanisms were observed, which were not intended for our study. In a pooled sample experiment, the contribution of these individual-specific reactions is diluted, while the general mechanisms active in all animals remain emphasized.

Finally, the molecular interaction networks of DEGs were investigated using the Cytoscape v. 3.1.0 software (http://www.cytoscape.org/) (Saito et al. 2012) as a complementary and more comprehensive approach for the identification of central hub genes. The networks of gene relationships were based on the Homo sapiens gene annotation databases. Results were visualized using the ClueGO v. 2.1.1. app. of the Cytoscape software to create clusters of functionally related genes, using the gene ontology databases: KEGG, Wiki Pathways, REACTOME, GO Immune System Process, GO Molecular Function, GO Cellular Component, and GO Biological Process.

The probability value was calculated using Enrichment/Depletion (a two-sided hypergeometric test) with Bonferroni correction. The networks of genes were created for the following situations: (1) all DEGs, (2) DEGs expressed only in the E group and were below the detection level for the C group, (3) DEGs expressed only in the C group and were below the detection level in the E group, (4) DEGs expressed only in the C or E group when combined in one network. The parameters “number of genes per node” and “minimum percentage of genes per node” were chosen empirically for different analyses.

Results

Animals and fatty acid profiles

There were no significant differences in growth performance and carcass characteristics between the pigs in the two treatment groups (data not shown). However, a tendency toward a lower daily gain and a lower backfat thickness in the experimental group was observed, but it was statistically insignificant. No interaction was observed between addition of PUFA and slaughter weight.

The analysis of the fatty acid profile showed a higher level of PUFAs in the liver of the experimental group (p-value = 0.0026), as a consequence of the increased content of PUFAs in the diet. Supplementation of the diet with PUFAs was associated with a higher content of most omega-3 fatty acids in the liver and unchanged or lower omega-6 fatty acid content, hence, pronounced reduction of the omega-6/omega-3 fatty acid ratio (Table 3). Particularly, contents of eicosapentaenoic acid (EPA) (p-value = 0.0006) and docosahexaenoic acid (DHA) (p-value = 0.0006) of the omega-3 family were higher, while that of arachidonic acid (AA)—omega-6 was lower in the E group in comparison to the C group (p-value = 0.0055). Interestingly, the mean value of the omega-6/omega-3 fatty acid ratio was approximately 30 in standard diet group (C group) and approximately 4 for the group with increased content of PUFAs.

Table 3 Fatty acid content (g/100 g FAME) in liver samples from two dietary groups. Data are presented as averages

We did not observe any statistically significant differences in regard to the content of saturated fatty acids (SFA). We observed a lower content of MUFA in the liver of the E group compared to the C group (p-value = 0.07).

Next-generation sequencing

Next-generation sequencing using RNA-Seq approach was performed to characterize the pig liver transcriptome. The number of reads per sample obtained from the NGS was 13,509,248 sequences using paired-end libraries (PE). Among these, 9,812,246 (72.63 %) reads were mapped in pairs, 381,668 (2.83 %) were mapped in broken pairs, and 3,315,334 (24.54 %) reads were not mapped. The average read was 151 base-pair lengths.

In regard to the fragment counting, the total number of fragments was 6,754,624, in which 1,848,501 (27.37 %) remained uncounted. Among 4,906,123 counted fragments (72.63 % of total fragments), 4,401,376 (65.16 %) were unique and 504,747 (7.47 %) were non-specific. An intact pair was counted as one, and broken pairs were ignored.

In total, 20,325 transcripts were annotated to genes revealed by the RNA-Seq technology for the pig liver samples (Table 4). Analysis identified 3565 DEGs between the two dietary groups, 1484 of which were up-regulated and 2081 of which were down-regulated in the E group compared to the C group. The differentially expressed genes identified in the both dietary groups included 2098 annotated transcripts. Results also indicated that expression of certain genes in the liver were unique to a particular dietary group. For instance, 809 genes were unique to the E group while the expression of these genes in the opposite group was below our detection level and 658 genes were unique to the C group.

Table 4 Summary of RNA sequencing analysis data of pig liver tissues fed a control (C) diet or the experimental (E) diet consisting of additional omega-6 and omega-3 fatty acids

The mean value of the absolute expression level was 67.95 for all reads in both dietary groups. However, high-throughput sequencing results also identified genes that were characterized by an extreme high absolute value of expression level (RPKM ≥1000). Among the 3565 DEGs, 30 genes in the E group and 41 genes in the C group had an extreme high absolute value of expression, 27 of which were common to both dietary groups.

Fold change (FC) for highly expressed genes (RPKM ≥1000) ranged mainly from 2.0 to 5.0. Two genes had higher FC—retinol saturase (RETSAT, FC = 9.1) and serine dehydratase (SDS, FC = 27.3). Among the up-regulated genes characterized by the highest value of absolute expression were as follows: APOA4, APOA5, HMGCS2, ACSL1, and CYP2C49, while highly expressed down-regulated genes were COX1, SULT2A1, HPD, FMO1, and CPS. These genes are mainly related to fatty acid activation, fatty acid β-oxidation, lipid binding, and lipid transport activity (up-regulated) or were enzymes involved in prostaglandin biosynthesis or peroxidase activity (down-regulated).

The mean FC value of all DEGs was 9.4. The next-generation sequencing analysis also revealed genes characterized by a high FC value caused by a low level of expression for genes within one dietary group.

Quantitative real-time PCR

Differentially expressed genes selected from RNA-Seq generated data were validated using quantitative real-time PCR. Almost all genes chosen for validation were regulated in the same direction in both approaches excluding EXTL1 and FASN. The differences in Log2-fold of change between both methods were observed as shown in Fig. 2a. The differences in mean relative expression between control and experimental groups were observed (Fig 2b) in agreement with RNA-Seq data using RPKM expression values for almost all genes excluding EXTL1 (opposite results) FASN and FADS2 (statistically insignificant). High variability between the particular individuals among groups was visible.

Fig. 2

Validation of differentially expressed genes (DEGs) by qRT-PCR. a comparison of the RNA-Seq (6 animals per group) and qRT-PCR (15 animals per group), Log2-fold of change between experimental and control groups, b individual variability of validated differentially expressed genes in qRT-PCR. For qRT-PCR, fold change was assessed on the basis of mean values of gene expression for 2 dietary groups with 15 individual samples per group. All Ct values were normalized to the topoisomerase (DNA) II beta (TOP2B) gene. Error bars represent standard error of variability among individual samples. *p-value <0,05; **p-value <0,01; ***p-value <0,001, ns not significant

Functional categorization of DEGs and pathway analysis

Fatty acid metabolism

The different ratio of omega-6/omega-3 fatty acids in the liver was reflected in altered fatty acid metabolism pathways revealed by the Cytoscape-ClueGo software. The visualization of potential networks for selected genes was consistent with previous observations of the most highly expressed genes and revealed the major biological mechanisms related to fatty acid metabolism: triglyceride metabolic process, lipid homeostasis, the PPAR signaling pathway, and fat digestion and absorption (Fig. 3 and Additional file 1).

Fig. 3

The gene network pathways related to fatty acid metabolism. Visualization generated by Cytoscape-ClueGo software after selection criteria: at least three genes per one pathway node and 4 % of the associated genes from all genes in one node

All DEG

To further investigate the biological importance of DEG, we performed gene category analysis in DAVID software. Among gene sets comprising all DEGs, 1623 records were assigned to gene identifiers recognized by the DAVID tool. Table 5 shows 12 significant functional annotation clusters (for more details see Additional file 2). The created clusters represented genes involved in cell membrane-related processes and signal transduction. The cluster with the highest enrichment score (8.6) included 363 genes related to the plasma membrane. The functional clusters that demonstrated the highest number of involved genes were transmembrane proteins in the cellular component category, signal peptide and glycoproteins in the biological processes category, and ion binding for the molecular function category. Taken together, major functions within this group appeared to be related to plasma membrane proteins, extracellular matrix, and cell transport and signaling. Functional classification in the DAVID tool adjusted DEG to terms in the KEGG database and showed the importance of 3 KEGG pathways—cytokine-cytokine receptor interactions (42 genes), focal adhesion (34 genes), and dilated cardiomyopathy (21 genes).

Table 5 Functional clusters created for all DEGs

Cytoscape bioinformatics analysis of potential interactions between proteins for all DEGs was performed on 729 genes (Additional file 3). The major biological functions revealed by Cytoscape-ClueGo software were as follows: positive regulation of the JAK-STAT cascade, muscle cell development, dilated cardiomyopathy, growth of symbiont in the host, positive regulation of leukocyte proliferation, regulation of intracellular protein transport, positive regulation of inflammatory response, response to molecules of bacterial origin, and cerebellum development (Fig. 4). Several of these pathways were grouped around biological themes, such as inflammatory/pathogen response (10 pathways), signaling pathways (4), brain/nerve development (3), carbohydrate/lipid metabolism (3), and muscle development and function (2). The JAK-STAT pathways (on the right, near the immune response pathways) and the pathogen response pathway (in the lower right) were grouped relatively close to the inflammatory response pathways (Fig. 4). Muscle development was located near the top of the figure, above and close to the cardio-related pathways. The main groups representing the brain development pathways are in the upper left, separated from the other pathways. The carbohydrate/lipid metabolism pathway is on the bottom left. The regulation of the intracellular protein transport pathway is located in the center. All signaling-related pathways were dispersed throughout the network. Overall, the processes, which were grouped together, were more related to each other than they were to the others. The distribution of pathways for all DEGs as specific clusters is presented in Fig. 5.

Fig. 4

The gene network and biological processes for all DEGs. Figure generated by Cytoscape-ClueGo. The red and green dots represent the up-regulated and down-regulated genes in the E group compared to the C group, respectively. The lines represent interactions among genes. The processes marked in grey are processes without direct connections to any of the other processes. The selection criteria for this network were to have at least two genes per node with a minimum of 5 % of the associated genes from all loaded genes in one node. The group percentage was 50 % with a kappa score threshold of 0.4

Fig. 5

The pathway distribution of the specific clusters for all DEGs. Analysis performed in Cytoscape-ClueGo. The selection criteria for this network were to have at least two genes per node with a minimum of 5 % of the associated genes from all loaded genes in one node. The group percentage was 50 % with a kappa score threshold of 0.4

Genes expressed only in the experimental group

The analysis in DAVID, which comprised set of genes specific to the E group, recognized and analyzed 359 gene IDs. The results for this subset of genes are presented in Table 6 (more details are listed in the Additional file 4). The cluster with the highest enrichment score was the glycoprotein/transmembrane cluster. The signal peptide and ion transport clusters belong to the functional mechanisms in this category. No pathways were revealed from the KEGG database for the set of genes expressed exclusively in the E dietary group.

Table 6 Functional clusters for genes expressed only in E diet group

Visualization of the molecular interaction networks in Cytoscape for DEGs detected only in the E group was performed for 59 genes (Additional file 5). Selected genes were involved mainly in the regulation of steroid hormone biosynthesis (Fig. 6). Other affected pathways involved muscle contraction, nervous system, signaling, immune response, and ion channels. The results also included a number of small networks that were separated from the main network.

Fig. 6

The gene network and biological processes for genes specific to experimental (E) group. Analysis performed in Cytoscape-ClueGo. The selection criteria for this network were to have at least two genes per node with 15 % of the associated genes from all loaded genes in one node. The kappa score threshold was 0.51

Genes expressed only in the control group

In total, 275 DEGs were converted to the DAVID-recognized gene IDs, but no functional cluster was assigned in this group (Additional file 6). The KEGG pathway analysis in DAVID software revealed two molecular interaction networks: cytokine-cytokine receptor interaction (12 genes) and type I diabetes mellitus (5 genes).

Cytoscape classification of DEGs specific for the C group was conducted on 31 genes (Additional file 7). Adjusted pathways included chemokine biosynthetic process, inflammatory response, and regulation of calcidiol 1-monooxygenase activity involved in vitamin D biosynthesis (Fig. 7). Besides the main network, other functions located outside of this network included potassium-transporting ATPase activity, cGMP biosynthetic processes, and syncytium formation.

Fig. 7

The genes network and biological processes for genes specific to control (C) group. Analysis performed in Cytoscape-ClueGo. The selection criteria for this network were to have at least two genes per node with 15 % of associated genes from all loaded genes in one node. The kappa score threshold was 0.51

Genes exclusively expressed in both groups

Cytoscape analysis preformed on DEGs identified exclusively in the E or C group but not in both of the groups showed that the most of the pathways were related to neurological system processes (5 pathways), immunity response (2), and cation transport (2) (Fig. 8). Two specific clusters were assigned in the interaction network processes: regulation of neurological process and regulation of alcohol biosynthetic process.

Fig. 8

The genes network for genes expressed above detection only in one of the groups merged in one network. Analysis performed in Cytoscape-ClueGo. The data sets of genes expressed only in one of groups were combined together in one analysis to create one figure. The selection criteria for this network were to have at least two genes per node with a minimum of 4 % of the associated genes from all loaded genes in one node. The group percentage was 50 % with a kappa score threshold of 0.4

Discussion

In our study, the RNA-Seq analysis showed modulatory effects of dietary omega-6 and omega-3 fatty acids on mRNA expression in pig liver. Liver transcriptomes of pigs fed a diet with increased LA and ALA content (the precursors for the synthesis of long-chain PUFA) were compared to those of pigs fed the standard diet. The scope of our study was to investigate the general biological mechanisms underlying the physiological effects of increased dietary intake of PUFA.

Increase of omega-3 PUFA contribution in fatty acid profile

As we expected, our data indicated specific differences between dietary groups in regard to the fatty acid profile. Results showed a higher contribution of PUFAs among all analyzed series of fatty acids in liver for pigs fed the LA and ALA-enriched diet compared to the standard diet, particularly for the omega-3 fatty acids EPA and DHA, but not omega-6. These results are consistent with those reported for previous gene expression studies in pigs (Corominas et al. 2013) and also rats (Hanke et al. 2013). This finding indicates higher absorption, transport and storage, or lower degradation of unsaturated fatty acids, which suggests inhibition of lipogenesis, but not for the omega-3 pathway. Substantial differences in the omega-6/omega-3 ratio in the liver, lowered ratio in the tissue after PUFA supplementation, outlined that the diet presumably affected body composition, which was the basis for further investigation on the transcriptome level. We can suppose that in the group characterized by increased content of LA and ALA in the diet, appeared specifically increased absorption of omega-3 fatty acids, not omega-6, on the base of fatty acid profile. Even though we consider that the conversion of ALA to EPA and DHA in vivo is less than 2 % (Brenna et al. 2009; Burdge and Calder 2005), the effect on omega-3 pathway was clear.

Enhancement of fatty acid β-oxidation in liver through induction of the PPAR signaling pathway

Among the most highly DEGs (FC ≥ 8), we identified genes involved in lipid metabolism, which play a key role in the synthesis of fatty acids. We observed lower expression of those genes in the liver of pigs with increased omega-3 fatty acids contribution in the liver, suggesting a major down-regulation of lipid metabolism and increase of β-oxidation. This is thought to relate to the omega-3 and omega-6 fatty acids (Jump 2011). One of the identified genes characterized by significantly reduced expression in our study was CYP7A1, which was shown to be important for the PPAR signaling pathway and is one of the target genes of LXR. The liver X receptor plays a major role in the PPAR signaling pathway and is a key regulator of fatty acid synthesis by regulation of lipogenic gene transcription and induction of SREBP-1c gene expression (Afman and Muller 2012). The effects of omega-3 fatty acids on changes in the expression of LXR’s target genes agreed with the findings of a recent study in pigs (Ramayo-Caldas et al. 2012).

Additionally, the expression of another transcription factor, PPARα, was affected. The peroxisome proliferator-activated receptor is a liver-specific gene and plays a central role in fatty acid β-oxidation in hepatocytes. It is thought to be the main indicator of the fatty acid oxidation-metabolic state (Gormaz et al. 2010). A higher content of omega-6 and omega-3 fatty acids in our study resulted in higher expression of the PPARα gene in the liver. Ramayo-Caldas et al. (2012) reported a similar effect on nuclear factors where the influence on the liver transcriptome of a PUFA-rich diet was compared to monounsaturated fatty acid (MUFA)- and saturated fatty acid (SFA)-enriched diet. Here, we pointed out how different contents of omega-6 and omega-3 PUFAs in the diet resulted in a different ratio of these two fatty acid families in the liver and how these alterations are reflected in the transcriptome. In our study, it is proposed that an omega-6 and omega-3 PUFA-enriched diet induced β-oxidation on the transcriptomic level in comparison to the influence of standard diet. It is generally accepted that omega-3 fatty acids induce hepatic fatty acid oxidation, possibly via PPARα activation, and suppress hepatic fatty acid synthesis through SREBP-1c (for reviews see El-Badry et al. (2007). However, some investigators have reported contrary results. Hanke et al. (2013) showed that omega-3-rich diets had no apparent effect on molecular markers of fatty acid oxidation on protein level. In particular, PPAR and ACOX were not affected nor were markers of fatty acid synthesis (SREBP-1c and ACC1) for livers of rats with diet-induced obesity.

Lastly, we found that up-regulated genes of the genes with the highest value of absolute expression registered by RNA-Seq analysis were involved in fatty acid activation, fatty acid β-oxidation, lipid binding, and lipid transport activity. To summarize, the direction of the expression level for major transcription factors responsible for fatty acid homeostasis, such as PPARα, LXR, or SREBP-1c, showed that omega-6 and omega-3-enriched diets influenced transcription of fatty acid β-oxidation genes and had significant impact on their molecular status what highlight their role as bioactive components of a diet. The higher PPARα expression in the liver suggested the biological modulatory effect of the LA- and ALA-enriched diet at the molecular level.

Increase of conversion intensity for omega-3 PUFA pathways in the liver

In our study, the FADS2 and FADS1 genes, which are responsible for coding enzymes delta-6-desaturase and delta-5-desaturase, respectively, where up-regulated in the RNA-Seq analysis (validated by qRT-PCR, statistically significant for FADS1). They are common for the conversion of fatty acids on three different pathways of omega fatty acid metabolism, being involved in introducing additional double bonds to dietary PUFAs, linoleic (omega-6) and α-linolenic acid (omega-3), and also to oleic acid (omega-9 MUFA) derived from de novo synthesis or the diet (Ratnayake and Galli (2009).

On the base of transcriptomic analysis, we cannot unambiguously confirm to which pathway delta-6-desaturase and delta-5-desaturase are directed. However, we can hypothesize a potential mechanism from the contributions of the derivatives of omega-6 and omega-3 fatty acid precursors in the fatty acid profile in the liver. Increased content of EPA and DHA (omega-3) and decreased AA (omega-6) suggest that ALA is converted into EPA and further to DHA, but conversion of LA to AA seems to be suppressed. Therefore, we suggest an increased intensity of conversion for omega-3 PUFA in the liver and suppose that FADS2 and FADS1 act mainly in the omega-3 metabolic pathway.

Furthermore, we observed a decreased expression of lipogenic gene coding for acetyl-CoA carboxylase-1 (ACACA). This suggests that the synthesis of short-chain fatty acids was decreased as a consequence of influence of the diet on these lipogenic genes. However, one of another lipogenic gene fatty acid synthase (FASN) was not confirmed by qRT-PCR. Conversely, enhanced activity of the conversion occurs to the long-chain fatty acids related to increased expression of ELOVL6 gene coding elongase (confirmed by qRT-PCR), which convert SFA-palmitic acid to stearic acid and is limiting enzyme for conversion between SFA and MUFA.

The production of acetyl-CoA derived from fatty acid degradation processes does not impact the Krebs cycle independently from the PUFAs content in the diet.

The main difference between the dietary groups in our studies was increased omega-3 intake from the diet and consequently higher conversion of omega-3 in the fatty acid pathway metabolism. Increased contents of omega-6 and omega-3 fatty acids in the diet resulted in altered metabolic pathways in the liver in a specific manner. The main biochemical pathways of unsaturated fatty acid metabolism seem to be directed to the omega-3 pathway, but not omega-6.

Further considerations in our future studies will concern identification of genes participating in the regulation of enzymes involved in lipid metabolism, which are common for different pathways of fatty acid metabolism.

Activation of pathways related to signal transduction, membrane composition, and mitigation of inflammation state

Analysis in DAVID demonstrated that most of the differentially expressed genes identified among two treatment groups affected by increased omega-3 fatty acid conversion were involved in signal transduction and plasma membrane-related proteins (Table 5). The biological function of signaling pathways is to forward information from the extracellular region via receptors on the plasma membrane and ion channels to support their transfer to the cytoplasm, what is supposed to be influenced by a low ratio of omega-6/omega-3 PUFA in our study and highlighted by the DAVID analysis. For example, extracellular fibronectins, one of the functional clusters revealed by DAVID, are ligands for the integrins located on the plasma membrane surface (Huveneers et al. 2008). An interaction network of reactions mediate between the transmission of various extracellular signals to the nucleus and is responsible for positive regulation of molecular function. All biological functions listed in the functional annotation clustering from DAVID are related to the processes taking place between cell membrane and the extracellular matrix such as ion transport or kinases activation. Polyunsaturated fatty acids are important components of the cell membranes and influence membrane fluidity as well as the action of membrane-bound enzymes and receptors. Thus, it is not surprising that oral administration of omega-6 and omega-3 fatty acids in pigs resulted in alterations emphasizing major functions of fatty acids in organisms—signal transduction and building of membrane, revealed in DAVID and Cytoscape-ClueGo analysis.

DAVID analysis concatenated with the KEGG database further support findings mentioned above and revealed focal adhesion molecules as a significant pathway, which is responsible for signal transduction, functionality of cell membranes, kinases activity, and ion transport. This pathway was related to the PUFA supplemented group, for which the omega-3 fatty acid content was increased.

Additionally, processes specific to pigs fed the supplemental diet (E group) were analyzed, and we observed results consistent with our findings mentioned above concerning the directions of alterations in lipid metabolism. We observed that the regulation of steroid hormone biosynthesis was the most affected pathway of the up-regulated genes, which were expressed only in the group characterized by increased intensity of conversions on omega-3 pathway. Steroid hormones are breakdown products of the EPA and DHA conversion pathways, and these metabolites were influenced by the administration of excess omega-3 fatty acid precursors in the diet.

Conversely, down-regulated genes, expressed only in the standard diet group, represent signaling molecules assigned as a cytokine-cytokine receptor interaction pathway identified by the KEGG database in the DAVID, (see also below) what reflects the anti-inflammatory properties of a high content of omega-3 fatty acids or low ratio of omega-6/omega-3 fatty acids in the liver. The expression of pro-inflammatory cytokines, which occurs in the control group of pigs under a normal physiological condition, is not observed in the PUFA supplemented group, what can hint at mitigation of the inflammation state. Our findings are consistent with studies on omega-3 fatty acid interventions against inflammatory status (Liu et al. 2015). Omega-3 fatty acids have been shown to inhibit expression of cytokine mRNAs, which are stimulated via omega-6 fatty acids (Kruger et al. 2010; Roy et al. 2007). Omega-3 fatty acids disturb the action of omega-6 fatty acid excess. Thus, our results suggest that the diet enriched with omega-3 fatty acids can modulate the immune response and could possibly improve the health of animals.

When we analyzed the processes specific to pigs fed the standard diet, we also observed perturbations of health-related pathways. Under a normal, physiological state, in the standard diet group of pigs, cytokine interaction was the most prominent pathway for genes expressed above detection level registered only in this group, while in the group with increased conversion of omega-3 fatty acids, expression of these genes were either not expressed or were down-regulated. Therefore, we can assume that increased conversion of omega-3 fatty acids revealed anti-inflammatory effect. This is in agreement with a recent study that investigated the effects of different omega-6/omega-3 ratios in pigs (Duan et al. 2014). Duan et al. (2014) revealed that inhibition of immune stimulation in pigs was characterized by a low ratio of omega-6/omega-3 fatty acids. Similarly, Papadopoulos et al. (2009) observed inflammation modulation in sows after administration of high doses of omega-3 fatty acids during the perinatal period. The up-regulation of TNFα in our study, which is an indicator of the inflammatory state, and additionally KEGG pathway and Cytoscape results further support this observation.

Additionally, calcidiol 1-monooxygenase activity, which is involved in vitamin D synthesis, was down-regulated. Vitamin D is known to promote vascular regeneration and possesses anti-proliferative and pro-apoptotic functions. Deficiency of vitamin D was associated with the inflammatory state (Wong et al. 2014). Presence of pro-inflammatory cytokines (such as TNFα and IL-6) might impair activation of vitamin D, which plays a pivotal role in the regulation of the prostaglandin pathway, thereby affecting both PGE2 synthesis by inhibition of COX2 expression (down-regulated in our studies), as well as its degradation by stimulation of 15-PGDH expression (Hummel et al. 2014). Thus, it could be proposed that a higher conversion of omega-3 fatty acids in the liver, caused by increased dietary intake of LA resulted in inhibition of inflammation in vascular tissue possibly via activation of the vitamin D system.

In summary, the presented work utilized high-throughput sequencing to identify the biological and metabolic processes of the pig liver transcriptome most affected by different omega-6 and omega-3 fatty acid content in the diet. A low ratio of omega-6/omega-3 fatty acids in the liver, caused by increased intake of omega-3 fatty acids from the diet, modulated gene expression and affected a set of genes related mainly to energy metabolism, signaling pathways, and the inflammatory response.