Skip to main content

Advertisement

Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency

Article metrics

Abstract

Background

Selection for feed efficiency is crucial for overall profitability and sustainability in dairy cattle production. Key regulator genes and genetic markers derived from co-expression networks underlying feed efficiency could be included in the genomic selection of the best cows. The present study identified co-expression networks associated with high and low feed efficiency and their regulator genes in Danish Holstein and Jersey cows.

RNA-sequencing data from Holstein and Jersey cows with high and low residual feed intake (RFI) and treated with two diets (low and high concentrate) were used. Approximately 26 million and 25 million pair reads were mapped to bovine reference genome for Jersey and Holstein breed, respectively. Subsequently, the gene count expressions data were analysed using a Weighted Gene Co-expression Network Analysis (WGCNA) approach. Functional enrichment analysis from Ingenuity® Pathway Analysis (IPA®), ClueGO application and STRING of these modules was performed to identify relevant biological pathways and regulatory genes.

Results

WGCNA identified two groups of co-expressed genes (modules) significantly associated with RFI and one module significantly associated with diet. In Holstein cows, the salmon module with module trait relationship (MTR) = 0.7 and the top upstream regulators ATP7B were involved in cholesterol biosynthesis, steroid biosynthesis, lipid biosynthesis and fatty acid metabolism. The magenta module has been significantly associated (MTR = 0.51) with the treatment diet involved in the triglyceride homeostasis. In Jersey cows, the lightsteelblue1 (MTR = − 0.57) module controlled by IFNG and IL10RA was involved in the positive regulation of interferon-gamma production, lymphocyte differentiation, natural killer cell-mediated cytotoxicity and primary immunodeficiency.

Conclusion

The present study provides new information on the biological functions in liver that are potentially involved in controlling feed efficiency. The hub genes and upstream regulators (ATP7b, IFNG and IL10RA) involved in these functions are potential candidate genes for the development of new biomarkers. However, the hub genes, upstream regulators and pathways involved in the co-expressed networks were different in both breeds. Hence, additional studies are required to investigate and confirm these findings prior to their use as candidate genes.

Background

Globally, food demand is increasing as a consequence of world population growth [1]. However, arable land to produce sufficient amounts of food is decreasing, and the carbon footprint is increasing [2]. Hence, solutions for efficient and environmentally friendly methods to produce food are urgently needed.

Feed efficiency (FE) in dairy cattle is the ability of a cow to convert the feed nutrient consumed into milk and milk by-products. Many approaches have been developed and adopted to select the most feed-efficient cows. Currently, residual feed intake (RFI) has been used to measure FE in dairy cows [3, 4]. Residual feed intake is the difference between the predicted and actual feed intake [5]. Regression models have been used to calculate the RFI value. Thus, animals with low RFI values are more efficient [6]. The genetic selection of animals with a low RFI will improve profitability [7], decrease greenhouse gasses emissions [8] and optimize the use of food resources. However, in the case of dairy cattle, the interpretation of RFI is not straightforward. Many other factors should be considered, as this selection might lead to a negative energy balance, cause health issues and affect the fertility of the cows [9, 10].

In Denmark, Holstein and Jersey are the most common dairy breeds used [11]. Comparatively, Holstein and Jersey cattle do not differ in terms of digestibility, energy efficiencies, and the ability to convert dietary protein to milk protein [12]. However, there are no gene expression profiling studies of these breeds. Hence, to understand the complex biological mechanisms in nutrient partitioning in dairy cattle, liver transcriptomics analysis may be useful to interpret and understand the pathways and functional elements of the genomes involved [13]. Transcriptomics is a form of high throughput analysis to quantify gene expression in a specific cell type or tissue [14]. Various studies have reported that mRNA levels of many genes are heritable, which affects genetic analysis [15,16,17]. Many studies based on transcriptomics (microarray and RNA-sequencing) have been conducted to study gene expression in feed efficiency [18,19,20]. Studies on differential gene expression have been well established to identify candidate genes for biomarker development [21]. There are limited studies related to gene expression for RFI traits in dairy cattle, particularly for Jersey and Holstein breeds. However, some studies have reported the gene expression associated with RFI in other breeds and species. For example, Lkhagvadorj et al. [22] found that the common energy consumption controlled by PPARA, PPARG and/or CREB is related to RFI in pigs. In beef cattle, Alexandre et al. [19] reported the alteration of lipid metabolism and an increase in the inflammatory response in animals with low feed efficiency. Paradis et al. [20] also reported a greater response to hepatic inflammation in heifers with high feed efficiency. In Nellore beef cattle, Tizioto et al. [23] identified the differentially expressed genes involved in oxidative stress. Hence, transcriptomics analysis might provide additional knowledge on the complex mechanisms that regulate nutrient intake.

Diet affects the energy metabolism and efficiency of dairy cows [24]. Some studies have investigated the correlation between FE and diet, focusing on the gene expression profiles of specific tissues. Dairy cows are typically fed high energy or high-concentrate feed to meet the high-energy demand during the lactation period. It has previously been shown that high energy feeding does not affect the fatty acid concentration but does affect the expression of genes such as ACACA, LPL and SCD in the lipid metabolism [25]. Thus, it is also interesting to investigate the effects of different levels of energy in feed using co-expression network approaches.

Previously, we performed differential gene expression analysis on RNA from the livers of Holstein and Jersey cows. We identified several differentially expressed genes between high and low RFI [26]. The differentially expressed genes were related to primary immunodeficiency, steroid hormone biosynthesis, retinol metabolism, starch and sucrose metabolism, ether lipid metabolism, arachidonic metabolism and cytochrome P450 in drug metabolism. These biological processes and pathways are important mechanisms that are associated with feed efficiency.

Therefore, it is important to thoroughly investigate the mechanisms controlling feed efficiency. Systems biology is the most promising approach to obtain a better understanding of complex traits, such as feed efficiency. In systems biology, many computational methods are based on network approaches. Co-expression network analysis has been successfully used to analyse complex traits and diseases in humans and animals [27,28,29,30]. Weighted Gene Co-expression Network Analysis (WGCNA) can be used to identify clusters (modules) of highly correlated genes [31]. WGCNA has been used to identify candidate genes that are associated with the FE. Alexandra et al. (2015) identified differentially co-expressed genes that are involved in lipid metabolism in RFI divergent Nellore cattle. Similarly, lipid metabolism-related processes were identified in low-RFI pigs [22].

In the present study, the WGCNA method was applied to RNA-Seq data from the livers of Holstein and Jersey cows to: i) identify groups of co-expressed genes and biological pathways associated with RFI; ii) identify the hub genes and upstream regulators in these modules that may be good candidate genes for feed efficiency-related traits; and iii) compare the mechanisms and processes involved in RFI between Holstein and Jersey cattle. To our knowledge, this study is the first to use weighted gene network approaches to examine the overall complex transcriptional regulation of feed efficiency (RFI) using RNA-Seq data in Danish Holstein and Jersey cows.

Materials and methods

Animal ethics statement

The experimental design and animals that were being used in this experiment were permitted by the Danish Animal Experimentation Inspectorate.

Experimental data

The experimental design and details of the experimental animals have been previously described in [26].

In brief, the dataset used in this experiment consists of 38 RNA-Seq expression profiles of liver bioposies from nine Holsteins and ten Jersey cows. In each breed group, cows were classified in high and low feed efficient and RNA samples were collected before and after treatment diet (low and high concentrate diet). The animals were assigned to the different diets after at least for 14–26 days adaptation period. All 38 RNA samples were paired-end sequenced using Illumina HiSeq 2500. The bioinformatics pipeline for RNA-Seq data processing is described in [26]. The expression quantification was performed using Ensembl Bovine annotation (release 82). The raw count data matrix used in this study is available in http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92398.

Weighted gene co-expression network analysis (WGCNA)

The Weighted Gene Co-expression Network Analysis (WGCNA) [31] R package was used to build co-expression networks and identify groups of highly co-expressed genes. Individual analyses were conducted on each breed group.

First, the low count genes and outliers were filtered by leaving only genes that had at least 1 count per million in 90% of the group. The remaining 11,153 genes in Holstein and 11,238 genes in Jersey were used for the analysis. The gene expression counts were normalized using the default procedure from the DESeq2 package version 1.12.0 [32] by correcting for the parity number to reduce potential effects from the parity number factor. The normalized data were subsequently log transformed as suggested in the WGCNA manual (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/). The final dataset was used in WGCNA to build an unsigned network. Pairwise Pearson’s correlations among all genes were calculated to create an adjacency matrix. A soft threshold power was set at β = 12 for Holstein and β = 10 for Jersey, correspondent to a scale-free topology index (R2) [33] of 0.9 for Holstein and 0.8 for Jersey. The adjacency matrix was used to calculate the Topological Overlap Measure (TOM). Modules of co-expressed genes were identified by using the dynamic tree cut algorithm [34]. Modules were arbitrarily labelled with different colours.

The module eigengenes were computed for each module using the first principal component to capture the variation in gene expression within each module. The eigengene sign was chosen to have a positive correlation with average module gene expression.

The correlation between module eigengene and RFI or treatment diet was evaluated to select modules that were associated with the respective traits (p-value < 0.05). In addition, FDR were computed using Benjamini–Hochberg (BH) method separately for each breed.

Gene significance (GS) was computed for each gene as the correlation between gene expression counts and FE. In addition, hub genes were identified, selecting genes with high module membership (MM > 0.8) in the modules of interest.

Functional enrichment analysis

The modules that are significantly associated with RFI and treatment diet traits were selected.

Functional enrichment analysis was performed in the selected modules to identify and interpret complex biological functions based on gene ontology terms for the biological processes, molecular functions and cellular components and based on the KEGG pathways annotation.

All the genes included in each module were used in the functional enrichment analysis with the Cytoscape 3.4.0 plug-in software, ClueGO v2.2.6 [35]. The significance value was set as p-value < 0.05 and the BH correction was used as the multiple test correction. The reference set used for this analysis included a total of 9064 genes. The list of genes in the module of interest was also analysed using the STRING v.10.0 [36] database and the Bos taurus annotation.

Ingenuity® Pathway Analysis (IPA®) was used to detect upstream regulators, diseases and functions in the selected modules. The upstream regulator analysis identifies the upstream regulators that better explain the change in gene expression. The analysis is based on the set of indirect relationships present in the IPA® database. The algorithm computes an overlap P-value by measuring enrichment of network-regulated genes to determine the most likely set of upstream regulators. Next, the algorithm computes the activation Z-score by identifying the match of up- and down-regulation annotated in Ingenuity knowledge base. The Z-score is then used to predict the activation state of the upstream regulators (either activated or inhibited).

A summary of the pipeline of the experimental workflow, bioinformatics and statistical analysis is presented in Fig. 1.

Fig. 1
figure1

Experimental design and co-expressed gene network analysis pipeline

Results

In the present study, WGCNA was used to identify RFI and diet-associated co-expression modules and their key functions. In total, 72 modules (Fig. 2) for Holstein cows and 59 modules (Fig. 3) for Jersey cows were identified. Subsequent the module detection, we have performed multiple testing corrections (Additional file 1: Tables S1 and S2 in each breed using BH method despite the norm that it is not carried out across gene network modules and traits. Unfortunately, after the multiple testing corrections, none of the top module is significant at adjusted p-value < 0.05 and therefore the results are to be validated in independent experiments with larger sample size, which is beyond the scope of this study. The results reported here are therefore are of exploratory and preliminary in nature. Therefore, modules with nominal p-value< 0.05 were used to be reported and discussed in the subsequent sections.

Fig. 2
figure2

Module trait relationship (p-value) for detected modules (y-axis) in relation with traits (x-axis) for Holstein cows. The module trait relationship were colored based on the correlation between the module and traits (red = strong positive correlation; green = strong negative correlation). X-axis legend: Diet = Treatment diet; RFI = Residual feed intake; Lact_no = Lactation number

Fig. 3
figure3

Module trait relationship (p-value) for detected modules (y-axis) in relation with traits (x-axis) for Jersey cows. The module trait relationship were colored based on the correlation between the module and traits (red = strong positive correlation; green = strong negative correlation). X-axis legend: Diet = Treatment diet; RFI = Residual feed intake; Lact_no = Lactation number

A total of 11 modules and four modules were significantly correlated with RFI for Holstein and Jersey cows, respectively. Additionally, 13 modules for Holstein and two modules for Jersey were significantly associated with treatment diet.

We assigned all the significant modules into the ClueGO application analysis to investigate the gene ontology (GO) and KEGG pathway-related functions with specific traits. The modules with the top significant module trait relationships (MTRs) were selected as the modules of interest in the present study. The modules lightsteelblue1 and violet in Jersey cows and the modules salmon and magenta in Holstein cows were selected for RFI and treatment diet, respectively.

Modules related to RFI and treatment diet in Holstein cows

In Holstein cows, among the 11 modules that were significantly (p-value< 0.05) related to the RFI, salmon module (203 genes with MTR RFI = 0.7) is the top significant module. For the diet trait, we identified the magenta module as the top significant module. The magenta module comprised 212 genes that contribute to the MTR Diet = 0.82.

In the top module (salmon), steroid biosynthesis was identified as the most enriched KEGG pathway (Fig. 4). This finding was also confirmed after analysing the genes in this module using www.string-db.org, and almost the same pathways and same patterns appeared in the output. Interestingly, most of the enriched pathways of co-expressed genes in Holstein cows were involved in steroid, lipid and cholesterol biosynthesis and metabolism (Fig. 4).

Fig. 4
figure4

Pie chart presenting an overview of the significant GO terms and KEGG pathways in the salmon module in Holstein cows

Additional file 1: Table S3 shows a summary of the functional groups with the number of genes involved in the GO terms and pathways. In total, 84 GO terms were significantly enriched (p-value< 0.05) after multiple testing corrections using BH. The GO-terms and KEGG pathways presented here are also almost the same as the output from the STRING 10 analysis (Additional file 1: Tables S5, S6 and S7).

The list of upstream regulators identified for the modules that are significantly associated with RFI and diet are presented in Additional file 1: Table S11. In the salmon module, ATP7B was predicted as activated, while POR and cholesterol were predicted as inhibited. In Additional file 1: Tables S13 and S14 shows the diseases and functions involved in salmon and magenta modules.

The module eigengene diagram for both of the salmon and magenta modules shows a higher average expression profile in high RFI samples (Fig. 5 a and b).

Fig. 5
figure5

a Module eigengene (y-axis) across samples (x-axis) from the salmon module (associated to RFI) (b) Module eigengene (y-axis) across samples (x-axis) from the magenta module (associated to treatment diet)

The list of genes with high (MM > 0.8) in the salmon module is presented in Table 1.

Table 1 List of the top hub genes generated from (MM > 0.8) in the salmon module in Holstein cows

Modules related to RFI and treatment diet in Jersey cows

Among the four modules significantly (p-value< 0.05) related to RFI in the Jersey group, the lightsteelblue1 module (72 genes) with a module trait relationship (MTR RFI = − 0.57) is the top significant (p-value< 0.05) module associated with RFI. In total, 44 GO terms were significantly enriched (p-value< 0.05) after multiple test correction using BH. For the diet trait, among the two significantly correlated modules, the violet module was the top significant (MTR Diet = − 0.47). However, this module has limited output from a functional enrichment analysis or no interesting biological information related to diet. Hence, the modules related to diet for the Jersey breed were not further discussed.

Figure 6 and Additional file 1: Table S4 shows the top summarized GO terms involved in the lightsteelblue1 module that is related to immune system functions. The first and the second GO terms, which are associated with the regulation of lymphocyte activation and positive regulation of leukocyte activation, involved almost the same genes as those that are involved in immune system functions. In detail, primary immunodeficiency has been identified (p-value< 0.05) as a significant KEGG pathway that involves four genes together with the positive regulation of leukocyte activated GO terms.

Fig. 6
figure6

Pie chart visualization of GO terms and KEGG pathways in the lightsteelblue1 module in Jersey cows

We identified IFNG (Interferon Gamma) as inhibited and IL10RA (Interleukin 10 Receptor Subunit Alpha), NKX2–3 (NK2 Homeobox 3) and dexamethasone were predicted as activated upstream regulators (Additional file 1: Table S12). In Additional file 1: Tables S14 and S16 shows the diseases and functions involved in lightsteelblue1 and violet modules.

Interestingly, all of these upstream regulators have functions related to the immune system. In addition, GO-terms and KEGG pathways from the STRING 10 analysis (Additional file 1: Tables S8, S9 and S10) also give almost the same output.

The module eigengene for the lightsteelblue1 module shows has an average expression profile that is lower in high RFI individuals (Fig. 7).

Fig. 7
figure7

Module eigengene (y-axis) across samples (x-axis) from the lightsteelblue1 module (associated to RFI)

The list of genes with high (MM > 0.8) in the lightsteelblue1 module is presented in Table 2.

Table 2 List of the top hub genes generated from (MM > 0.8) in the lightsteelblue1 module in Jersey cows

Discussion

WGCNA identified groups of co-expressed genes that are expected to perform the same biological functions and affect RFI. From the MTR, we tested the modules that were significantly correlated to the focus traits (RFI and diet). However, only the most significant module had any interesting biological meaning associated with the traits (one module in each breed). Hence, only the most biologically meaningful modules were further analysed and discussed.

For Holstein cows, we identified pathways and upstream regulators related to steroid biosynthesis, lipid metabolism, cholesterol metabolism and production in salmon module. In particular, we identified the activation of cholesterol and lipid synthesis in high RFI cows. There was a tendency for these three mechanisms to be activated in the datasets, which is consistent with the idea that high synthesis of fat is correlated with the loss of energy used in milk production in dairy cows, resulting in less feed efficient animals [37]. This finding is also consistent with previous studies that associated high fat deposition with high RFI animals [6, 38]. The magenta module was significantly associated with diet and involved the energy consumption and regulation of glucose.

For Jersey cows, the lightsteelblue1 module was enriched for immune system-related functions. Interestingly, the upstream regulators for the genes in the lightsteelblue1 module (IFNG and IL10RA) were also related to the immune system. In particular, the immune system in high RFI group was activated. Thus, the activation of the immune system leads to low feed efficiency, which is consistent with previous studies [19, 39].

These findings are supported by evidence from the co-expression network analysis of both breeds.

Co-expressed networks in Holstein cows

The functional enrichment analysis determined that the module identified in Holstein cows was involved in cholesterol biosynthesis, steroid biosynthesis, lipid biosynthesis and fatty acid metabolism.

From the most significant pathways, cholesterol biosynthesis has previously been discussed, as its related genes are important in the RFI. The cholesterol biosynthetic pathway is responsible for the variability of cholesterol levels in cells [40]. This module was also enriched for lipid biosynthesis. Interestingly, the levels of cholesterol and lipids have previously been positively associated with RFI in beef cattle [41].

Many genes in this modules have previously been associated with feed efficiency, [39]. For example, Acetyl-CoA carboxylase alpha (ACACA), Acetyl-CoA Acetyltransferase 2 (ACAT2), and fatty acid synthase (FASN) genes in the modules are key genes in cholesterol biosynthesis, organic hydroxy compound metabolism, collagen fibril organization, steroid biosynthesis, astral microtubule organization, protein oligomerization and oxidoreductase activity, acting on the CH-CH group of donors and NAD or NADP as an acceptor. ACACA and FASN were found to be differentially expressed and co-expressed in other feed efficiency-related studies [22, 39, 42]. The main function of FASN is to catalyse the synthesis of palmitate from acetyl-CoA and malonyl-CoA, in the presence of NADPH, into long-chain saturated fatty acids. Hence, these genes have a tendency to affect the feed efficiency in Holstein cows. In addition, many studies have discussed the involvement of several genes included in the modules that we identified in the present study (CYP7A1, ACACA, FASN) [39, 43]. The presence of ACAT2 is also interesting because the product of this gene is involved in lipid metabolism [44].

Other feed efficiency studies, for example, in pigs, have previously observed that lipogenesis and steroidogenesis in liver tissue are closely related to feed efficiency [22, 45], confirming previous observations in the differential expression analysis of this dataset [26].

In Holstein cows, we identified ATP7B as a top upstream regulator for the salmon module. This protein uses energy in the molecule adenosine triphosphate (ATP), which is responsible for the transport of metals into and out of cells using the energy stored in the molecule adenosine triphosphate (ATP). ATP7B appears to be activated in high RFI (low FE). Hoogeveen et al. (1995) [46] stated that the deficiency of copper in rats would increase the utilization of fat in rats. Hence, this finding suggests a relationship when ATP7B is activated, which potentially reflects the deposition of fats. Consistent with the present study, the high RFI cow shows the activation of ATP7B. This upstream regulator shows a relationship with regulating the fat consumption. Although it is not straightforward, the presence of the gene reflects the consumption of fat and indirectly affects the fat composition [47].

In the present study, cholesterol synthesis was activated in the IPA upstream regulator analysis. Furthermore, the activation of lipid metabolism in the disease function analysis supports the evidence from the GO term and pathway analyses. As lipid and cholesterol metabolism, and fat synthesis in particular, are activated in the high RFI group, we can assume that the high RFI group is inefficient in converting fat to energy. Hence, animals with high RFI (low FE) have high levels of cholesterol and fat in the body [48]. This finding is also consistent with Arthur et al. [49], who reported the positive relationship between RFI and average back fat in beef carcasses.

Interestingly, when fed a high or low concentrate diet, triglyceride homeostasis was the top GO biological process, which might be the result of the high energy or low energy diet. A previous study reported that controlled diet (with fructose and glucose) significantly affects the triglyceride levels [50].

Generally, based on the results obtained from the functional enrichment analysis for the Holstein breed, the most important GO terms, KEGG pathways and upstream regulators involved were related to steroid biosynthesis, cholesterol biosynthesis, lipid biosynthesis and triglyceride homeostasis. These findings show that the feed efficiency in Holstein cows is strictly associated with the regulation of energy via lipid and cholesterol metabolism.

Co-expressed networks in Jersey cows

The most significant pathways in Jersey cows were positive regulation of interferon-gamma production, lymphocyte differentiation, side of membrane, natural killer cell-mediated cytotoxicity and primary immunodeficiency. Interestingly, these most summarized pathways were related to the immune system. From the IPA upstream regulator and diseases function analysis, the immune system related functions were activated in the high RFI group.

Several studies also suggested that the involvement of the immune system would affect the feed efficiency [51, 52]. For example, [19, 27] discussed important findings but in different species and breeds. Kristina et al. [39] discovered an increase in the inflammatory response of the progeny of low RFI sires, which is consistent with the results of the present study. The type of diet might also affect the immune response. For example, Ametaj et al. [53] reported that the feeding of high concentrate feeds affects several inflammatory responses in feedlot steers. However, in the present study, no significant effect from the different type of concentrate diet in Jersey cows was observed. This finding might reflect the different populations and different breeds, as dairy cattle convert their nutrients into different products with respect to beef cattle [54]. Although, many other studies relate their findings with the importance of the immune system in RFI and feed efficiency, few studies have been conducted in dairy cattle [19, 20, 39].

Furthermore, these significant GO terms and pathways were also supported by the findings from upstream regulator analysis through IPA®. The top upstream regulator in Jersey cows is Interferon Gamma (IFNG), which has an interesting relationship to interactions among nutrition, metabolism, and the immune system [55]. This gene encodes a soluble cytokine that is a member of the type II interferon class. IFNG was predicted to be inhibited in high RFI Jersey cows. This protein is secreted from cells of both the innate and adaptive immune systems. IFNG is important in the system because it directly inhibits viral replication. The down-regulation of this cytokine in the high RFI group in Jersey cows might affect the feed efficiency. Thus, IFNG plays an important role in regulating immune systems in animals. Another interesting upstream regulator in Jersey cows is IL10RA (Interleukin 10 Receptor Subunit Alpha), which was predicted to be activated. IL10RA is a receptor with anti-inflammatory properties [56]. The activation of this gene might result in inhibition of the synthesis of pro-inflammatory cytokines. Reynolds et al. [57] reported that IL10RA was differentially expressed in rumen papillae of divergent average daily gain steers and these authors showed a negative association between the inflammatory response and feed efficiency. Thus, the activation of IL10RA in the high RFI group would reflect the inflammatory response in Jersey cows.

We further speculate that, based on the results obtained in the Jersey breed, the most important GO terms, KEGG pathways and upstream regulators were related to the immune system. Jersey cows have many co-expressed genes that relate to the immune system to regulate feed utilization. It is likely that in Jersey cows, immunity plays a key role in substituting feed nutrient into milk and milk components. The immune response plays an important role in energy balance during milk production in dairy cows.

Comparison of RFI associated modules between Holstein and Jersey cows

In the datasets analysed in the present study, the most significant module associated with RFI differed between the Holstein and Jersey breeds. Furthermore, these modules were enriched for different sets of biological processes. This evidence suggests that the Holstein cow system is more reactive towards steroid biosynthesis, while Jersey cows have more reactions in their immune systems. Several studies have reported the importance of the lipid and cholesterol metabolism and immune system related functions in feed efficiency traits in farm animals, likely reflecting the complex role of the liver in regulating the nutrient uptake [58].

The hub genes of the modules identified in the present study represent potential candidate genes for RFI. These findings might provide additional information and new insights into the biological processes that are associated with RFI in these two main dairy breeds. Thus, we speculated that in this study population, the liver transcriptomics profiles of the two main dairy breeds are involved in two different biological processes. However, a comparative feed efficiency study reported similar results in terms of digestibility and ratios of milk to body weight and feed intake between Holstein and Jersey cows [12]. The sample size of the present study did not enable confirmation of whether the identified biological processes are breed specific. To confirm this notion, the set of genes should be validated in other cows using qPCR to confirm whether the expression patterns conform to different RFI-diet groups, which is out of the scope of the present study. In addition, a common limitation of static gene co-expression studies is the impossibility to decide if the identified modules are causing variation in the trait analysed or if their co-expression is a consequence of the variation observed for trait. Consequently, in this study we never talk about causality. Further study such as eQTL mapping (data integration between transcriptomics and genomics) could help in understanding better causality between gene expression and trait variation.

Conclusion

In conclusion, the co-expression network analysis revealed important genes and pathways in the liver that are involved in feed efficiency (RFI). In Holstein cows, the overall results showed that genes and upstream regulators such as ATP7b in RFI-associated modules that were co-expressed were primarily related to steroid and lipid biosynthesis. The results show that high RFI Holstein cows have a high lipid and cholesterol metabolism. The co-expressed genes associated with treatment diet were involved in triglyceride homeostasis. We observed different patterns of co-expressed genes involved in Jersey cows for which most of the co-expressed genes associated with RFI were related to the immune system in the most significant module. The upstream regulators IFNG and ILR10 that were predicted to be inhibited and activated, respectively, were closely associated with the immune system in Jersey cows. A high RFI Jersey cow tends to have a higher response to inflammation. The information of the functional enrichment from the analysis of co-expressed genes provides a better understanding of the mechanisms controlling RFI in Holstein and Jersey cows. Thus, the present study paves the way for the development of biomarkers for feed efficiency in dairy cattle.

Abbreviations

BH:

Benjamini-Hochberg

cDNA:

Complementary DNA

DM:

Dry Matter

FE:

Feed Efficiency

GO:

Gene Ontology

GS:

Gene Significance

HC:

High Concentrate

IPA®:

Ingenuity® Pathway Analysis

KEGG:

Kyoto Encyclopedia of Genes and Genomes

LC:

Low Concentrate

MM:

Module Membership

mRNA:

Messenger RNA

MTR:

Module Trait Relationship

RFI:

Residual Feed Intake

RNA-seq:

RNA-sequencing

TOM:

Topological Overlap Measure

WGCNA:

Weighted Gene CO-expression Network Analysis

References

  1. 1.

    Carolan M. Food security and policy; 2017.

  2. 2.

    Food, Nations AOotU. Food wastage footprint: impacts on natural resources: summary report. Rome: FAO; 2013.

  3. 3.

    Connor EE, Hutchison JL, Norman HD. Estimating feed efficiency of lactating dairy cattle using residual feed intake. Feed Efficiency in the Beef Industry; 2012. p. 159–73.

  4. 4.

    Varga G, Dechow C. Can we use residual feed intake to enhance dairy production efficiency? In: Proceedings of the 22nd Tri-State Dairy Nutrition Conference, Fort Wayne, Indiana, USA, 23–24 April 2013: 2013: Michigan State University; 2013. p. 131–40.

  5. 5.

    Koch RM, Swiger LA, Chambers D, Gregory KE. Efficiency of feed use in beef cattle. J Anim Sci. 1963;22(2):486–94.

  6. 6.

    Herd R, Arthur P. Physiological basis for residual feed intake. J Anim Sci. 2009;87(14_suppl):E64–71.

  7. 7.

    Williams Y, Pryce J, Grainger C, Wales W, Linden N, Porker M, Hayes B. Variation in residual feed intake in Holstein-Friesian dairy heifers in southern Australia. J Dairy Sci. 2011;94(9):4715–25.

  8. 8.

    Arthur JP, Herd R. Residual feed intake in beef cattle. Rev Bras Zootec. 2008;37(SPE):269–79.

  9. 9.

    Veerkamp R, Koenen E, De Jong G. Genetic correlations among body condition score, yield, and fertility in first-parity cows estimated by random regression models. J Dairy Sci. 2001;84(10):2327–35.

  10. 10.

    Llewellyn S, Fitzpatrick R, Kenny D, Murphy J, Scaramuzzi R, Wathes D. Effect of negative energy balance on the insulin-like growth factor system in pre-recruitment ovarian follicles of post partum dairy cows. Reproduction. 2007;133(3):627–39.

  11. 11.

    RYK- Ydelseskontrol, https://www.landbrugsinfo.dk/Kvaeg/RYK/Sider/Avereagemilkyieldisnowcloseto10000kgprcowandyears.aspx, Accessed date: 12 June 2017.

  12. 12.

    Blake RW, Custodio AA, Howard WH. Comparative feed efficiency of Holstein and Jersey Cows1. J Dairy Sci. 1986;69(5):1302–8.

  13. 13.

    Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.

  14. 14.

    Suravajhala P, Kogelman LJ, Kadarmideen HN. Multi-omic data integration and analysis using systems genomics approaches: methods and applications in animal production, health and welfare. Genet Sel Evol. 2016;48(1):38.

  15. 15.

    Li J, Burmeister M. Genetical genomics: combining genetics with gene expression analysis. Hum Mol Genet. 2005;14(suppl 2):R163–9.

  16. 16.

    Brem RB, Yvert G, Clinton R, Kruglyak L. Genetic dissection of transcriptional regulation in budding yeast. Science. 2002;296(5568):752–5.

  17. 17.

    Schadt EE, Molony C, Chudin E, Hao K, Yang X, Lum PY, Kasarskis A, Zhang B, Wang S, Suver C. Mapping the genetic architecture of gene expression in human liver. PLoS Biol. 2008;6(5):e107.

  18. 18.

    Kadarmideen HN. Systems biology in animal production and health, vol. 2. Switzerland: Springer; 2016.

  19. 19.

    Alexandre PA, Kogelman LJ, Santana MH, Passarelli D, Pulz LH, Fantinato-Neto P, Silva PL, Leme PR, Strefezzi RF, Coutinho LL. Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle. BMC Genomics. 2015;16(1):1073.

  20. 20.

    Paradis F, Yue S, Grant J, Stothard P, Basarab J, Fitzsimmons C. Transcriptomic analysis by RNA sequencing reveals that hepatic interferon-induced genes may be associated with feed efficiency in beef heifers. J Anim Sci. 2015;93(7):3331–41.

  21. 21.

    Zhao X-M, Qin G. Identifying biomarkers with differential analysis. In: Bioinformatics for diagnosis, prognosis and treatment of complex diseases. The Netherlands: Springer; 2013. p. 17–31.

  22. 22.

    Lkhagvadorj S, Qu L, Cai W, Couture OP, Barb CR, Hausman GJ, Nettleton D, Anderson LL, Dekkers JC, Tuggle CK. Gene expression profiling of the short-term adaptive response to acute caloric restriction in liver and adipose tissues of pigs differing in feed efficiency. Am J Phys Regul Integr Comp Phys. 2010;298(2):R494–507.

  23. 23.

    Tizioto PC, Coutinho LL, Oliveira PS, Cesar AS, Diniz WJ, Lima AO, Rocha MI, Decker JE, Schnabel RD, Mourão GB. Gene expression differences in longissimus muscle of Nelore steers genetically divergent for residual feed intake. Sci Rep. 2016;6:39493.

  24. 24.

    Armentano L, Weigel K. Considerations for improving feed efficiency in dairy cattle, vol. 2013; 2013. p. 37.

  25. 25.

    Tao H, Chang G, Xu T, Zhao H, Zhang K, Shen X. Feeding a high concentrate diet down-regulates expression of ACACA, LPL and SCD and modifies milk composition in lactating goats. PLoS One. 2015;10(6):e0130525.

  26. 26.

    Salleh M, Mazzoni G, Höglund J, Olijhoek D, Lund P, Løvendahl P, Kadarmideen H. RNA-Seq transcriptomics and pathway analyses reveal potential regulatory genes and molecular mechanisms in high-and low-residual feed intake in Nordic dairy cattle. BMC Genomics. 2017;18(1):258.

  27. 27.

    Kogelman LJ, Cirera S, Zhernakova DV, Fredholm M, Franke L, Kadarmideen HN. Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA sequencing in a porcine model. BMC Med Genet. 2014;7(1):57.

  28. 28.

    Kadarmideen HN, Watson-Haigh NS, Andronicos NM. Systems biology of ovine intestinal parasite resistance: disease gene modules and biomarkers. Mol BioSyst. 2011;7(1):235–46.

  29. 29.

    Kogelman LJ, Byrne K, Vuocolo T, Watson-Haigh NS, Kadarmideen HN, Kijas JW, Oddy HV, Gardner GE, Gondro C, Tellam RL. Genetic architecture of gene expression in ovine skeletal muscle. BMC Genomics. 2011;12(1):607.

  30. 30.

    Cogill SB, Wang L. Co-expression network analysis of human lncRNAs and cancer genes. Cancer Inform. 2014;13(Suppl. 5):49.

  31. 31.

    Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics. 2008;9(1):559.

  32. 32.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

  33. 33.

    Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4(1):1128.

  34. 34.

    Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the dynamic tree cut package for R. Bioinformatics. 2008;24(5):719–20.

  35. 35.

    Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman W-H, Pagès F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3.

  36. 36.

    Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2014;43(D1):D447–52.

  37. 37.

    Bauman DE, Currie WB. Partitioning of nutrients during pregnancy and lactation: a review of mechanisms involving homeostasis and homeorhesis. J Dairy Sci. 1980;63(9):1514–29.

  38. 38.

    Xi Y, Wu F, Zhao D, Yang Z, Li L, Han Z, Wang G. Biological mechanisms related to differences in residual feed intake in dairy cows. Animal. 2016;10(08):1311–8.

  39. 39.

    Weber KL, Welly BT, Van Eenennaam AL, Young AE, Porto-Neto LR, Reverter A, Rincon G. Identification of gene networks for residual feed intake in angus cattle using genomic prediction and RNA-seq. PLoS One. 2016;11(3):e0152274.

  40. 40.

    Parhami F, Mody N, Gharavi N, Ballard AJ, Tintut Y, Demer LL. Role of the cholesterol biosynthetic pathway in osteoblastic differentiation of marrow stromal cells. J Bone Miner Res. 2002;17(11):1997–2003.

  41. 41.

    Karisa B, Moore S, Plastow G. Analysis of biological networks and biological pathways associated with residual feed intake in beef cattle. Anim Sci J. 2014;85(4):374–87.

  42. 42.

    Graugnard DE, Piantoni P, Bionaz M, Berger LL, Faulkner DB, Loor JJ. Adipogenic and energy metabolism gene networks in longissimus lumborum during rapid post-weaning growth in Angus and Angus× Simmental cattle fed high-starch or low-starch diets. BMC Genomics. 2009;10(1):142.

  43. 43.

    Horodyska J, Oster M, Wimmers K, Mullen A, Lawlor P, Hamill R. P3024 transcriptome analysis of longissimus thoracis et lumborum from pigs divergent in residual feed intake. J Anim Sci. 2016;94(supplement4):63–4.

  44. 44.

    Nguyen TM, Sawyer JK, Kelley KL, Davis MA, Rudel LL. Cholesterol esterification by ACAT2 is essential for efficient intestinal cholesterol absorption: evidence from thoracic lymph duct cannulation. J Lipid Res. 2012;53(1):95–104.

  45. 45.

    Zhao Y, Hou Y, Liu F, Liu A, Jing L, Zhao C, Luan Y, Miao Y, Zhao S, Li X. Transcriptome analysis reveals that vitamin a metabolism in the liver affects feed efficiency in pigs. G3: Genes| Genomes| Genetics. 2016;6(11):3615–24.

  46. 46.

    Hoogeveen RC, Reaves SK, Lei KY. Copper deficiency increases hepatic apolipoprotein AI synthesis and secretion but does not alter hepatic total cellular apolipoprotein AI mRNA abundance in rats. J Nutr. 1995;125(12):2935.

  47. 47.

    Huster D, Purnat TD, Burkhead JL, Ralle M, Fiehn O, Stuckert F, Olson NE, Teupser D, Lutsenko S. High copper selectively alters lipid metabolism and cell cycle machinery in the mouse model of Wilson disease. J Biol Chem. 2007;282(11):8343–55.

  48. 48.

    Richardson E, Herd R, Archer J, Arthur P. Metabolic differences in Angus steers divergently selected for residual feed intake. Anim Prod Sci. 2004;44(5):441–52.

  49. 49.

    Arthur P, Archer J, Johnston D, Herd R, Richardson E, Parnell P. Genetic and phenotypic variance and covariance components for feed intake, feed efficiency, and other postweaning traits in Angus cattle. J Anim Sci. 2001;79(11):2805–11.

  50. 50.

    Schaefer EJ, Gleason JA, Dansinger ML. Dietary fructose and glucose differentially affect lipid and glucose homeostasis. J Nutr. 2009;139(6):1257S–62S.

  51. 51.

    Van Eerden E, Van Den Brand H, Parmentier H, De Jong M. Phenotypic selection for residual feed intake and its effect on humoral immune responses in growing layer hens. Poult Sci. 2004;83(9):1602–9.

  52. 52.

    Mpetile Z. Effect of divergent selection for residual feed intake on immune system of Yorkshire pigs; 2014.

  53. 53.

    Ametaj B, Koenig K, Dunn S, Yang W, Zebeli Q, Beauchemin K. Backgrounding and finishing diets are associated with inflammatory responses in feedlot steers. J Anim Sci. 2009;87(4):1314–20.

  54. 54.

    Roberts A, Funston R, Mulliniks T, Petersen M, MacNeil M. Feed efficiency-how should it be used for the cow herd? 2011.

  55. 55.

    Schroder K, Hertzog PJ, Ravasi T, Hume DA. Interferon-γ: an overview of signals, mechanisms and functions. J Leukoc Biol. 2004;75(2):163–89.

  56. 56.

    Gondret F, Vincent A, Houée-Bigot M, Siegel A, Lagarrigue S, Louveau I, Causeur D. Molecular alterations induced by a high-fat high-fiber diet in porcine adipose tissues: variations according to the anatomical fat location. BMC Genomics. 2016;17(1):120.

  57. 57.

    Reynolds J, Foote A, Freetly H, Oliver W, Lindholm-Perry A. Relationships between inflammation-and immunity-related transcript abundance in the rumen and jejunum of beef steers with divergent average daily gain. Anim Genet. 2017;48(4):447–9.

  58. 58.

    Jungermann K, Katz N. Functional specialization of different hepatocyte populations. Physiol Rev. 1989;69(3):708–64.

Download references

Acknowledgements

We gratefully acknowledge collaborator scientists, specifically Peter Lund, Johanna Höglund and Dana Olijhoek, for contributions to the present study. We also thank the field staff and laboratory technicians who assisted with the experiments in the present study.

Funding

This research was funded by the “Feed Utilization in Nordic Cattle” from Danish Milk Levy Foundation, Skejby Denmark. SMS received a stipend for PhD studies from the Ministry of Higher Education of Malaysia and Universiti Putra Malaysia and a tuition fee waiver from the University of Copenhagen. The funding bodies played no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

The data discussed in this publication were deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE92398 at: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92398

Author information

HNK was the overall project leader who conceived and conducted this selective transcriptomics profiling study on cows with high/low feed efficiency, and supervised SMS in the laboratory work and bioinformatics/systems biology analyses. PL provided feed efficiency measurements for the cattle in the experiment. SMS processed liver tissue for the RNA isolation and quality control of the RNA samples prior to RNA Sequencing. SMS analysed the data and drafted the original manuscript with assistance from GM. GM substantially contributed to the bioinformatics analyses. All authors wrote, read, and approved the final version of the manuscript.

Correspondence to H. N. Kadarmideen.

Ethics declarations

Ethics approval and consent to participate

The experimental design and animals that were being used in this experiment were permitted by the Danish Animal Experimentation Inspectorate.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1 and S2. Multiple testing corrections for the module trait relationship using Benjamini-Hochberg (BH) method. Table S3 and S4. ClueGO analysis output. Table S5-S10. STRING 10 analysis for salmon module in Holstein and lightsteelblue1 module in Jersey. Table S11 and S12. Upstream regulators from IPA® analysis. Table S13-S15. Disease and functions for most significant modules. (DOCX 65 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Salleh, S.M., Mazzoni, G., Løvendahl, P. et al. Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency. BMC Bioinformatics 19, 513 (2018) doi:10.1186/s12859-018-2553-z

Download citation

Keywords

  • RNA-seq
  • Feed efficiency
  • Residual feed intake
  • Co-expressed genes
  • Hub genes
  • Pathways
  • Holstein
  • Jersey
  • Dairy cattle