Skip to main content

Hepatic transcriptome perturbations in dairy cows fed different forage resources

Abstract

Background

Forage plays critical roles in milk performance of dairy. However, domestic high-quality forage such as alfalfa hay is far from being sufficient in China. Thus, more than 1 million tons of alfalfa hay were imported in China annually in recent years. At the same time, more than 10 million tons of corn stover are generated annually in China. Thus, taking full advantage of corn stover to meet the demand of forage and reduce dependence on imported alfalfa hay has been a strategic policy for the Chinese dairy industry. Changes in liver metabolism under different forage resources are not well known. Thus, the objective of the present study was to investigate the effect of different forage resources on liver metabolism using RNAseq and bioinformatics analyses.

Results

The results of this study showed that the cows fed a diet with corn stover (CS) as the main forage had lower milk yield, DMI, milk protein content and yield, milk fat yield, and lactose yield than cows fed a mixed forage (MF) diet (P <  0.01). KEGG analysis for differently expressed genes (DEG) in liver (81 up-regulated and 423 down-DEG, Padj ≤0.05) showed that pathways associated with glycan biosynthesis and metabolism and amino acid metabolism was inhibited by the CS diet. In addition, results from DAVID and ClueGO indicated that biological processes related to cell-cell adhesion, multicellular organism growth, and amino acid and protein metabolism also were downregulated by feeding CS. Co-expression network analysis indicated that FAM210A, SLC26A6, FBXW5, EIF6, ZSCAN10, FPGS, and ARMCX2 played critical roles in the network. Bioinformatics analysis showed that genes within the co-expression network were enriched to “pyruvate metabolic process”, “complement activation, classical pathway”, and “retrograde transport, endosome to Golgi”.

Conclusions

Results of the present study indicated that feeding a low-quality forage diet inhibits important biological functions of the liver at least in part due to a reduction in DMI. In addition, the results of the present study provide an insight into the metabolic response in the liver to different-quality forage resources. As such, the data can help develop favorable strategies to improve the utilization of corn stover in China.

Background

Forage is the largest component of the diet of lactating cows and could affect dry matter intake (DMI) [1] and consequently milk performance. Nutrient content such as crude protein (CP), neutral detergent fiber (NDF), and non-fibrous carbohydrate (NFC) of different forages differ greatly. For instance, alfalfa hay, a well-known high quality forage, contains higher CP, rumen degradable protein (RDP), and rumen undegradable protein (RUP) content than corn stover and Chinese wild rye grass [2]. Increasing RUP by 1% can improve milk production by 1 kg [3]. In addition, cows fed high proportions of alfalfa hay have higher milk protein production by increasing microbial protein yield, which may be attributed to the increased supply of rumen-available energy [2].

While high-quality forage such as alfalfa hay is still a bottleneck for the development of the dairy industry in China, it has become one of the largest dairy producers in the world [4]. In 2019, more than 1.2 million tons of alfalfa hay were imported to cover the shortage of high-quality forage in China [5]. At the same time, it is estimated that more than 10 million tons of corn stover are generated annually in China [6]. Thus, taking full advantage of crop residues such as corn stover to meet the demand of forage and reduce dependence on imported alfalfa hay has been a strategic policy for the Chinese dairy industry [2].

In the last 10 years, a large number of studies to evaluate the nutritive value of corn stover or Chinese wild rye grass have been conducted [2, 4, 7, 8]. For instance, Zhu et al. (2013) investigated the effect of different forage sources on lactation performance, microbial protein (MCP) synthesis, and N utilization efficiency in early lactation dairy cows [2]. Through studying metabolites from four biofluids (rumen fluid, milk, serum, and urine), Sun et al. (2015) elucidated the metabolic mechanisms of milk production affected by forage quality [4]. Zhang et al. (2014) evaluated the effects of diets with three different quality forage sources (alfalfa hay, L. chinensis and cornstalk) on the rumen microbiota of dairy cows [7].

In ruminants, liver contributes to more than 80% of the glucose produced via gluconeogenesis [9, 10]. In addition, liver is a critical hub for numerous physiological processes including lipid metabolism, amino acid metabolism, detoxification, and immune defense [11, 12]. Overall function and metabolism of the liver are sensitive to the plane of nutrition of the cows. For instance, Shahzad et al. (2014) demonstrate that the liver of cows fed a diet to meet 80% of estimated requirements had greater lipid and amino acid catabolic capacity and a more pronounced cellular inflammatory and endoplasmic reticulum stress response, while the liver of cows fed to meet or exceed requirements had a larger cell proliferation and cell-to-cell communication and greater activation of pathways/functions related to triacylglycerol synthesis [13]. Previous studies have been mainly focused on the effect of different forage resources on lactation performance and rumen fermentation, but simultaneous changes in liver metabolism under different forage resources are not well known. Thus, the objective of the present study was to investigate the effect of different forage resources on liver metabolism using RNAseq and bioinformatics analyses.

Results

Milk performance of cows fed different forage resources

As shown in Table 1, milk yield (30.5 vs. 23.1 kg/d, P <  0.01) and efficiency (1.47 vs. 1.32%, P <  0.01) was lower with CS than MF. In addition, DMI (21.4 vs. 17.4 kg/d, P <  0.01), milk protein content and yield (3.66 vs. 3.32%, P <  0.01;1.11 vs. 0.77 kg/d, P <  0.01), milk fat yield (1.34 vs. 1.02 kg/d, P <  0.01), and lactose yield (1.47 vs. 1.13 kg/d, P <  0.01) were all decreased by CS compared with MF.

Table 1 Milk yield and composition of lactating cows fed diets based on different forage sources

Differently expressed genes (DEG) and functional analysis

A total of 8582 unigenes were detected in the liver of dairy cows and 504 DEG (81 up- regulated and 423 down-DEG, Padj ≤0.05) were identified between cows consumed CS and MF (Additional File 1). Functional analysis for the DEG was performed using DIA, DAVID, and ClueGO.

The whole DIA output is available in Additional File 1. As shown in Fig. 1 where the perturbation in CS cows vs. MF cows on the main categories of the KEGG pathways in liver is summarized, all categories and subcategories were inhibited to different extents. For instance, “Metabolism” followed by “Genetic Information Processing” and “Environmental Information Process” were the most impacted. Within the most impacted category of “Metabolism” (Fig. 1), the subcategory “Glycan Biosynthesis and Metabolism” was the most impacted and was overall inhibited. Among the top 20 impacted pathways in liver tissue of CS compared with MF cows uncovered by the DIA, most of the pathways were inhibited (Fig. 2). Few pathways such as “Sulfur relay system”, “Vitamin B6 metabolism”, and “Glycosaminoglycan biosynthesis-keratan sulfate” were highly activated, and “Glycosaminoglycan biosynthesis-ganglio series” and “alpha-Linolenic acid metabolism” were modestly activated. Furthermore, among the top 20 most impacted pathways, approximately 25% were related to “Glycan Biosynthesis and Metabolism” with the pathway of “Glycosphingolipid biosynthesis – globo series” being the most impacted (Fig. 2).

Fig. 1
figure1

Summary of the main categories and sub-categories of KEGG pathways as results of the transcriptomic effect on liver tissue of corn stover (CS) compared to mixed forage (MF) as analyzed by the Dynamic Impact Approach. On the right are the bar denoting the overall impact (in blue) and the shade denoting the effect on the pathway (from green – inhibited – to red – activated). Darker the color larger the activation (if red) or inhibition (if green) of the pathway

Fig. 2
figure2

The 20 most impacted pathways in liver tissue of corn stover (CS) compared to mixed forage (MF) uncovered by the Dynamic Impact Approach. On the right are the bar denoting the overall impact (in blue) and the shade denoting the effect on the pathway (from green – inhibited – to red – activated). Darker the color larger the activation (if red) or inhibition (if green) of the pathway

Results of DAVID analysis are shown in Fig. 3 where KEGG and GO Biological Process (GO_BP) analysis were conducted. The GO_BP analysis revealed that 14 and 6 different terms (P ≤ 0.05) were enriched by downregulated DEG and upregulated DEG respectively. For the KEGG analysis, there were 8 terms enriched among DEG in total, with 5 terms enriched with downregulated DEG (P ≤ 0.05).

Fig. 3
figure3

Significantly enriched Gene Ontology Biological Process and KEGG pathways revealed by DAVID analysis of the transcripts up- (in red shade in the figure) or down- (in blue shade in the figure) regulated in liver tissue of corn stover (CS) compared to mixed forage (MF) cows. In vertical axis is the terms, in horizonal axis is the transformed FDR (−log10PValue)

The GO_BP analysis was also performed using ClueGO (Fig. 4). The results show that downregulated DEG were enriched to “pyruvate metabolic process”, “positive regulation of proteasomal protein catabolic”, “amide biosynthetic process”, and “regulation of multicellular organism growth”, while the upregulated DEG were enriched to “myeloid cell development”, “Schwann cell development”, and “negative regulation of small GTPase mediated signal transduction” (P ≤ 0.05).

Fig. 4
figure4

Functional annotation of DEGs using ClueGO. a: Enriched by downregulated DEGs; b: Enriched by upregulated DEGs. Each node is a Gene Ontology (GO) Biological Process term. The size of the nodes reflects the statistical significance of each term. Larger the node size, smaller the P-value. Different node colors represent different functional groups. The name of each group is given by the most significant term of the group. The nodes are grouped by similarity of their associated genes

Co-expression Network and Functional analysis.

Co-expression network analysis provides insights into the patterns of transcriptome organization and can reveal common biological functions among network genes [14]. The co-expression network analysis of this study was conducted using DEG with correlation > 0.9 and Padj < 0.01 (Additional file 4). The entire co-expression network is shown in Fig. 5, and the annotation information of the genes is available in Fig. 6. As shown in Fig. 5, the co-expression network revealed 7 genes (FAM210A, SLC26A6, FBXW5, EIF6, ZSCAN10, FPGS, ARMCX2) with higher degree and betweenness centrality (ranking in top 7, Additional file 5) than others, indicating a more critical role played by them in the network.

Fig. 5
figure5

Co-expression networks constructed using differently expressed genes (DEG) with absolute correlation ≥0.9 and adjusted p-value ≤0.01 by Cytoscape. The color of the nodes represents the fold change of the gene expressed in mixed forage (MF) compared to corn stover (CS). Upregulated genes are in red color, downregulated genes in blue color. Deeper the color, higher the fold changes. The size of the nodes represents the combined ranking of the degree and betweenness of the nodes (genes) in the network. Larger the size, higher the ranking. Color of the edges represent the correlation between the genes. Positive correlation is in red, negative correlation is in blue. The width of the edge represents significance of the correlation between the two genes. Larger the width, smaller the Padj

Fig. 6
figure6

Gene Ontology Biological Process (GO_BP) annotation for the whole co-expression network. The size of the nodes reflects the statistical significance of each term. Larger the node size, smaller the P-value. Different node colors represent different functional groups. The name of each group is given by the most significant term of the group. The nodes are grouped by similarity of their associated genes

Annotation information analysis for the genes within the co-expression network was performed using ClueGO and is shown in Fig. 6. Genes within the whole network were significantly enriched in “complement activation, classical pathway”, “retrograde transport, endosome to Golgi”, “positive regulation of proteasomal ubiquitin-dependent protein catabolic process”, “microtubule bundle formation”, “negative regulation of supramolecular fiber organization”, “viral genome replication”, “protein localization to microtubule cytoskeleton”, “ribonucleoprotein complex localization”, and “pyruvate metabolic process” (Fig. 6).

Discussion

Liver plays a central role in supporting the anabolic capacity of the mammary gland. Net hepatic glucose production (3.1 kg/d) of mid-to-late lactating cows is able to meet glucose required for milk lactose synthesis and maintenance [15, 16]. In addition, liver plays dominant roles in determining the ultimate quantity and pattern of metabolites available for milk synthesis [16]. Metabolic function and, thus, energy metabolism of liver responds to a variety of environmental stimuli including fasting or level of feed intake [17], diet composition and productive (physiological) state [15]. Although a number of studies have been conducted to assess effects of low-quality forage resources on lactation performance and rumen fermentation [2, 4, 7, 8], there are limited data on the response by important organs such as the liver. Thus, we used transcriptomics and bioinformatics in an effort to better capture genome-wide transcriptional responses of dairy liver to feeding low-quality forage (CS) versus high-quality forage (MF).

Feeding CS reduces Milk performance

Consistent with a previous study where milk yield of cows fed more alfalfa than those fed corn stover (P = 0.07) decreased [2], in this study milk yield was lower with CS than MF (Table 1). In addition, milk protein content and yield, milk fat yield, and lactose yield were all decreased by CS compared with MF (Table 1). Clearly, a large portion of the decreased milk performance in this study was mainly attributed to the lower DMI of CS compared with MF cows [18]. The study of Zhu et al. (2013) showed that corn stover compared with alfalfa led to lower OM degradability in the rumen (53.2 vs. 47.8%, P = 0.01) [2], suggesting longer retention time of undegraded fiber. Thus, the lower DMI of CS vs. MF cows in this study was likely caused by excess bulk in the rumen.

Pathways in liver were extensively inhibited in CS cows vs. MF cows

In this study, all categories and subcategories of the KEGG pathways in liver were overall inhibited to different extents in CS vs. MF cows (Fig. 1). Furthermore, among the top 20 impacted pathways, in liver tissue of CS compared with MF cows uncovered by the DIA, most of the pathways were inhibited (Fig. 2). Data for inhibited pathways indicated an overall downregulated metabolism in liver of CS compared with MF cows, which agrees with results of Sun et al. (2015) in which ruminal fluid and serum metabolite concentrations decreased with a low-forage compared with high-forage diet [4]. Thus, together the data imply a decreased overall metabolism level when low-quality forage is fed.

Low-quality forage inhibited glycan biosynthesis and metabolism

As shown in Fig. 1, the subcategory “Glycan Biosynthesis and Metabolism” was the most impacted and was overall inhibited. Furthermore, among the top 20 most impacted pathways, approximately 25% were related to “Glycan Biosynthesis and Metabolism” with the pathway of “Glycosphingolipid biosynthesis – globo series” being the most impacted (Fig. 2). Glycans are simple or complex polymers composed of monosaccharides [19], and mediate a wide variety of biological processes including cell growth and differentiation, cell−cell communication, immune response, pathogen interaction, and intracellular signaling events [20]. At a molecular level, glycans are often the first points of contact between cells, and they function by facilitating a variety of interactions both in cis (on the same cell) and in trans (on different cells) [21]. Thus, the high perturbation of glycan biosynthesis and metabolism in this study suggests a potential effect of low-quality forage on hepatocyte communication or growth and differentiation, which was also validated by the results of DAVID and ClueGO where the biological process of “cell-cell adhesion” and “positive regulation of multicellular organism growth” were significantly enriched among the downregulated DEG (Fig. 3 and Fig. 4).

Among the top 20 impacted pathways, “Glycosphingolipid biosynthesis – globo series” and “Glycosphingolipid biosynthesis – lacto and neolacto series” were highly inhibited in CS vs. MF cows (Fig. 2). In addition, “Glycosphingolipid biosynthesis – ganglio series” was also highly impacted, but the change in direction of the DEG involved in the pathway was not consistent, which is embodied in the modest direction of the impact (Fig. 2). However, it was evident that glycosphingolipid biosynthesis metabolism was overall inhibited by CS vs. MF in this study. Glycosphingolipids (GSLs) comprise a heterogeneous group of membrane lipids formed by a ceramide backbone covalently linked to a glycan moiety [22], and are classified based on their carbohydrate structure into six major series in vertebrates including gangliosides, lacto-, neolacto-, muco-, isoglobo-, and globo-series GSL [23]. D’Angelo et al. (2013) compiled published data indicating that GSL could modulate various aspects of the biology of the cell including apoptosis, cell proliferation, endocytosis, intracellular transport, cell migration and senescence, and inflammation [22]. Zhang et al. (2004) concluded that specific changes in composition and metabolism of GSL occur during cell proliferation, cell cycle phases, brain development, differentiation, and neoplasia in various cell types [24]. In addition, GSL form “microdomains” or “rafts” within the cell membrane, which move within the fluid bilayer as platforms for the attachment of proteins during signal transduction and cell adhesion [24]. Thus, in this study, the inhibited glycosphingolipid biosynthesis metabolism seems to offer further proof that the communication or growth and differentiation of hepatocytes was potentially inhibited by the low-quality forage diet. The significance of the perturbation at a deeper level could not be ascertained by the results of the present study.

Inconsistent with the above 4 pathways, “Glycosaminoglycan biosynthesis - keratan sulfate” was highly activated in CS cows vs. MF cows (Fig. 2). Keratan sulfate (KS) is one of the glycosaminoglycans (GAG), occurring as keratan sulfate proteoglycans on the cell surface and in the extracellular matrix [25]. Pomin (2015) concluded that GAG displays anti-inflammatory functions by activating leukocyte rolling along the endothelial surface of inflamed sites and also regulating chemokine action on leukocyte guidance, migration and activation [26]. The study of Vailati-Riboni et al. (2016) in transition cows demonstrated that feeding at 125% of nutrient requirements activated hepatic GAG synthesis pathways in under-conditioned cows, while it inhibited it in optimally-conditioned cows [27]. Thus, it was suggested that overfeeding of fatter cows may decrease the synthesis of anti-inflammatory compounds and consequently induce some detrimental effects [27]. Taken together, previous and present data suggest activation of “Glycosaminoglycan biosynthesis - keratan sulfate” in response to feeding CS as an anti-inflammatory response. This idea was also validated by DAVID analysis where “complement activation” and “Complement and coagulation cascades” were significantly enriched with up-DEG (Fig. 3). However, the underlying mechanisms could not be ascertained from results of the present study.

The biosynthesis of KS is often markedly altered in response to metabolic, pathologic, or developmental changes in tissues [28]. Davies et al. (1999) suggested that the expression of keratan sulfate is down-regulated in migrating corneal endothelial cells, while abundance on the cell surface returns when cells cease migration [29]. Thus, this suggests that KS has an anti-migration character. However, the anti-adhesive properties of KS were previously reviewed by Caterson and Melrose (2018) and Funderburgh (2000) [28, 30]. Thus, the exact function of KS as it relates to cell-cell adhesion in hepatocytes is difficult to ascertain with the available data. In the present study, the paradoxical effect of CS vs. MF on cell adhesion was also highlighted by results of DAVID, where both “cell-cell adhesive” and “negative regulation of cell-matrix adhesion” were significantly enriched by the down- and up-regulated DEG, while “cell adhesive” was significantly enriched among the up-regulated DEG (Fig. 3).

Low-quality forage inhibits amino acid metabolism

Metabolism of amino acids was overall inhibited by low-quality forage (Fig. 1). Among the top 20 impacted pathways, “Arginine biosynthesis”, “Selenoamino acid metabolism”, “beta-Alanine metabolism”, and “Tryptophan metabolism” were all inhibited (Fig. 2). Similar results were also revealed by ClueGO where “amide biosynthesis process” and “positive regulation of proteasomal protein catabolic” were significantly enriched by downregulated DEG (Fig. 4). Sun et al. (2015) studied metabolite profiles from four biofluids (rumen fluid, milk, serum, and urine) of cows fed different forage resources using metabolomics, with 55, 8, 28, and 31 significantly different metabolites identified in the rumen fluid, milk, serum, and urine, respectively [4]. These metabolites were involved in glycine, serine, and threonine metabolism; tyrosine metabolism; and phenylalanine metabolism [4]. Sun et al. (2016) in a subsequent urine metabolomics analysis demonstrated that Tyr metabolism and Phe, Tyr and Try biosynthesis pathways had the most variation when corn stover replaced alfalfa hay [31]. The study of Wang et al. (2018) showed that cows fed CS had lower absorbable Leu in the duodenum, which suggested this diet led to shortage of microbial Leu [8]. Sun et al. (2015) demonstrated that, under different quality forage resources, the concentrations of Phe and Tyr in rumen fluid exhibited lower fold-change values (0.54 and 1.19, respectively) than those in the serum (1.01 and 1.34, respectively), which implied that Phe and Tyr may be utilized more in the liver of cows fed high-quality forage than compared with low-quality forage [4]. Thus, we speculate that the inhibition of amino acid metabolism in CS vs. MF cows in this study was suggestive of an inhibited amino acid utilization in liver in the cows fed low-quality forage diet.

Co-expression network analysis

In the co-expression network, degree represents the number of connections of a node in a network and betweenness centrality is the number of times that a path passes through the node, which represents the influence this node exerts over other nodes and their potential interactions in the network [32]. Thus, both degree and betweenness centrality are measures of the function of a node in network connectivity [33]. As shown in Fig. 5, the co-expression network revealed 7 genes (FAM210A, SLC26A6, FBXW5, EIF6, ZSCAN10, FPGS, ARMCX2) with higher degree and betweenness centrality (ranking in top 7, Additional file 5) than others, indicating a more critical role played by them in the network.

Among the 7 genes, FAM210A (a mitochondrial gene) which had the highest degree has a crucial role in regulating bone structure and function [34]. SLC26A6 belongs to the solute carrier 26 family, and encodes a protein involved in transporting chloride, oxalate, sulfate and bicarbonate [35,36,37,38,39]. Thus, the inhibited expression of SLC26A6 indicated a decreased transporting ability of chloride, oxalate, sulfate and bicarbonate in liver of CS cows (Fig. 5). FBXW5 is a the TSC2 binding receptor of CUL4 E3 ligase complex [40]. Hu et al. (2008) demonstrated that FBW5 (FBXW5) promotes ubiquitination of tumor suppressor TSC2 by DDB1-CUL4-ROC1 ligase, and depletion of FBW5 stabilizes TSC2 [41]. Ha et al. (2014) demonstrated that intracellular accumulation of TSC2 inhibits the activity of mTOR and increase autophagy [40]. Thus, in the present study, the downregulated FBXW5 seems to imply an intracellular accumulation of TSC2 and consequently an increase in autophagy in liver of CS vs. MF cows. EIF6 is eukaryotic translation initiation factor. Depletion of eIF6 (using specific siRNA-mediated knockdown) in Mz-ChA-2 and TFK-1 cell lines inhibit cell proliferation and induced apoptosis [42], while EIF6 over-expression increases the motility and invasiveness of cancer cells [43]. However, in this study, “positive regulation of apoptotic process” was significantly enriched by downregulated DEG implying that apoptosis was not induced by low-quality forage diet (Fig. 4). In addition, “protein folding” was significantly enriched by downregulated DEG (Fig. 4) and the pathway “Aminoacyl-tRNA biosynthesis” was highly inhibited by CS vs. MF. Thus, combined with the inhibited EIF6, the results of the present study suggested an inhibited protein synthesis in liver of cows fed low-quality forage.

FPGS is a gene encoding the folylpolyglutamate synthetase enzyme. This enzyme has a central role in establishing and maintaining both cytosolic and mitochondrial folylpolyglutamate concentrations and, thus, is essential for folate homeostasis and the survival of proliferating cells [44, 45]. Folate plays an essential role in nucleotide biosynthesis and biological methylation reactions as a methyl donor [46, 47]. Consistent with the downregulation of FPGS (Fig. 5), in this study, “Folate biosynthesis” was also inhibited by CS vs. MF (Fig. 2). Taken together, the results of the present study suggested a decrease in folate homeostasis in cows fed low-quality forage. Thus, the downregulation of “nucleotide-excision repair” may be a downstream cascade reaction due to decreased folate homeostasis (Fig. 3). In addition, the downregulated ZSCAN10 and ARMCX2 with high degree and betweenness centrality were also unraveled by co-expression network analysis (Fig. 5), but the significance of the perturbation was unclear.

Annotation information analysis for the genes within the co-expression network was performed using ClueGO and is shown in Fig. 6. A potential explanation for the perturbation of all these biological processes is beyond the scope of the present study. However, it is noteworthy that almost all these biological processes are energy-requiring except for “pyruvate metabolic process”, which is highly associated with energy metabolism. Furthermore, “pyruvate metabolic process” was also significantly enriched by the whole downregulated DEG in CS vs. MF (Fig. 4), indicating that energy metabolism in liver was inhibited by low-quality forage. Taken together, the inhibited energy metabolism unraveled by this study was suggestive of a central role in the whole metabolic perturbation in liver. The DMI of CS cows was indeed 4 kg/d lower than the MF cows, which may account at least in part for the decreased energy metabolism in liver in CS vs. MF.

Conclusion

As in previous studies, feeding a low-quality forage reduces production performance, but also leads to marked alterations in the hepatic transcriptome. Among the unique biological pathways identified through bioinformatics analysis, glycan biosynthesis and metabolism and amino acid metabolism were highly inhibited when the low-quality forage diet was fed. Biological processes related to cell-cell adhesion, multicellular organism growth, and amino acid and protein metabolism also were downregulated. Co-expression network analysis indicated that the downregulated genes related to autophagy and translation played critical roles in the network. In addition, pyruvate metabolic process and other energy-requiring biological processes were enriched in the co-expression network. Collectively, results indicated that, compared to high-quality forage diet, low-quality forage could inhibit several basic cellular functions of the liver. Furthermore, the results of the present study provide an insight into the metabolic response in the liver to different-quality forage resources. As such, the data can help develop favorable strategies to improve the utilization of corn stover in China.

Methods

Experiment animals and management

The field experiment of this study was performed in Beijing ZhongDi Dairy Farm (Beijing, China) and the cows used in this study were all from this farm. Thirty-two healthy lactating Holstein cows (body weight, 550 ± 10 kg; days in milk, 55 ± 15; daily milk yield, 31 ± 2.30 kg, primiparous) were selected and divided into two groups based on average daily milk yield, body weight and days in milk. The two groups were randomly assigned to two diets: (i) mixed forage diet (MF), and (ii) corn stover diet (CS). Each treatment included 16 cows (n = 16). The two diets were formulated to meet their nutrient requirement (net energy lactation) according to NRC (2001) [3]. The ingredients and the chemical composition of the two experimental diets are shown in Table 2 and Table 3. Forage-to-concentrate ratio (F:C) of the two diets were all 64:36. The diets were mixed daily and fed ad libitum as total mixed ration. The diet was supplied thrice per day at 07:00, 14:00, and 20:00 h in an equal amount that allowed for 5–10% orts. Cows were milked thrice daily at 07:00, 14:00, and 20:00 h and had ad libitum access to water. The total duration of the experiment lasted for 14 weeks, including an acclimatization period of two weeks. All the two group cows were fed MF diet in acclimatization period and randomly allocated to MF or CS in trail period.

Table 2 Ingredients composition of experimental diets
Table 3 Chemical composition of experimental diets

Sample collection

Dry matter intake of each cow was recorded using the RIC system (Hokofarm Group, Netherlands). The offered total mixed rations were sampled twice per week, and the samples were pooled for each week and analyzed for DM, CP, NDF, ADF, EE, Starch, Ca, and P content as previously described [48]. Five cows were randomly selected from each treatment for liver tissue collection. The biopsy was conducted at the end of the 14 weeks as described by Gao et al. (2019) and Bu et al. (2017) [49, 50]. Tissue samples were washed with PBS buffer prepared using RNAase-free water, and then stored in liquid nitrogen immediately until RNA extraction. Health was monitored post-surgery by recording rectal temperature, milk yield, and feed intake daily for 7 days. Surgical clips were removed 7 days post-biopsy and the cows were placed back to their original barns in the farm to continue rear.

RNA extraction and sequencing

Total RNA was extracted with TRIzol reagent (Life technologies, US, Cat#74106) according to the manufacturer’s protocol. Integrity and concentration of total RNA were then assessed using a 2100 Bioanalyzer (Agilent Technologies, US) with the RNA 6000 Nano Kit (Agilent Technologies, US). Complementary DNA (cDNA) was synthesised and used to construct a library with the NEBNext Ultra RNA Library Prep Kit (NEB, E7530). The libraries were sequenced on the Illumina HiSeq2000 platform via 2 × 50-bp paired-end sequencing at BGI Tech Solutions Co., Ltd. (Shenzhen, China). The RNA-Seq, libraries were sequenced at BGI Tech Solutions Co., Ltd. (Shenzhen, China).

Quality analysis, mapping, and Transcriptome assembly

The reads containing adapter, poly-N and low-quality reads in the raw data were removed to obtain the clean reads. All the downstream analyses were based on the clean data with high quality. Reference bovine genome and gene model annotation files were downloaded from (ftp://ftp.ensembl.org/pub/release-89/fasta/bos_taurus/dna/) and (ftp://ftp.ensembl.org/pub/release-89/gtf/bos_taurus) respectively. Index of the reference genome was built using Bowtie2 v2.2.8. HISAT2 (v2.0.4) was used to align paired-end clean reads to the reference genome [51]. HISAT2 was run with ‘--rna-strandness RF’, other parameters were set as default. The StringTie (v1.3.1) was used to assemble the mapped reads of each sample in a reference-based approach [52].

Statistical analysis

Data analysis of milk performance and milk composition was performed using the GLM model in SAS (8.2; SAS Institute Inc., Cary, NC), with treatment used as fixed effect, cow used as random effect, and day used as repeated effect.

Analysis of differently expressed genes (DEG)

The mapped reads count table of each gene and each sample was used for the analysis of DEG with DESeq2 package (1.26.0) in R (3.6.1) as the standard workflow instructed [53]. The complete dataset of DEG was in Additional file 1.

GO and KEGG enrichment analysis

Dynamic Impact Approach (DIA) was used to perform the analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [54]. The dataset including Entrez Gene ID, FDR (adjusted P-value, Padj), expression ratio, and P-value was uploaded and Entrez Gene ID were used as background. Adjusted P-value < 0.05 were used as cut-off. The output of DIA was exported in Additional file 2.

The enrichment analysis of various database including KEGG pathways, Gene Ontology Biological process (GO_BP), Cellular components (GO_CC) and Molecular function (GO_MF) was run by Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 [55]. For this analysis, all the annotated transcripts that were detected (Entrez Gene ID) were used as background and three datasets were analyzed: 1) up-regulated differently expressed genes (DEG, log2 fold change > 0, Padj < 0.05) by CS vs MF; 2) down-regulated DEG (log2 fold change < 0, Padj < 0.05) by CS vs MF; and 3) both up- and down-regulated DEG (Padj < 0.05) of CS vs MF. Results were downloaded using both Functional Annotation Chart and Functional Annotation Clustering (Additional file 3). The enrichment analysis of Gene Ontology Biological process was also conducted and visualized using ClueGO (2.5.5) a plug-in of Cytoscape (3.7.2) [56].

Co-expression network construction and functional annotation

Read counts of each gene and each sample were normalized by median normalization using the EBSeq (3.10) R package. The correlation and correlation significance of every pair DEG (Padj < 0.01) was calculated using logarithmic matrix of read counts with the R package “psych” (1.8.12). Pearson was chosen as the correlation method. Then, the table containing the statistically significant correlations across the whole data set for every pair of DEGs was generated (Additional file 4). R package “igraph” (1.2.4.1) was used to calculate network statistics (degree and betweenness centrality for each node) as described by Contreras-Lopez et al. (2018) [33]. The dataset is available in Additional file 5. The dataset including correlation of the DEGs, the degree and betweenness centrality of the nodes, and the log2 fold change of the DEGs were uploaded to Cytoscape. The gene symbols were set as the node’s identifiers. The correlation, correlation significance, degree and betweenness, and log2 fold change were mapped to the edge color, edge with, node size, and node fill color respectively. The functional annotation analysis of the network was performed using the ClueGO application. The network and functional annotation results of the network were shown in Fig. 5 and Fig. 6.

Availability of data and materials

The datasets (additional files) supporting the conclusions of this article are available in the figshare, the hyperlink to datasets is https://doi.org/10.6084/m9.figshare.12082941.v2. All sequencing data (fastq files) generated in the present study are available in the NCBI Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/sra/) with accession number SRP256053.

Abbreviations

ADF:

Acid detergent fiber

ARMCX2:

Armadillo repeat containing X-linked 2

CUL4:

Cullin 4

CP:

Crude protein

CS:

Corn stover

DAVID:

Database for Annotation, Visualization and Integrated Discovery

DDB1:

Damage specific DNA binding protein 1

DEG:

Differently expressed genes

DIA:

Dynamic Impact Approach

DM:

Dry matter

DMI:

Dry matter intake

EE:

Ether Extract

EIF6:

Eukaryotic translation initiation factor 6

FAM210A:

Family with sequence similarity 210 member A

FBXW5:

F-box and WD repeat domain containing 5

FDR:

False discovery rates

FPGS:

Folylpolyglutamate synthase

GAG:

Glycosaminoglycans

GO_BP:

Gene Ontology biological process

GO_CC:

Gene Ontology cellular components

GO_MF:

Gene Ontology molecular function

GSLs:

Glycosphingolipids

KEGG:

Kyoto Encyclopedia of Genes and Genomes

KS:

Keratan sulfate

MCP:

Microbial protein

MF:

Mixed forage

Mz-ChA-2:

Human gallbladder cancer

NDF:

Neutral detergent fiber

NFC:

Non-fibrous carbohydrates

RDP:

Rumen degradable protein

ROC1:

Rotamase CYP 1

RUP:

Rumen undegradable proteins

SLC26A6:

Solute carrier family 26 member 6

TFK1:

Human extrahepatic bile duct carcinoma cell line

TSC2:

TSC complex subunit 2

ZSCAN10:

Zinc finger and SCAN domain containing 10

References

  1. 1.

    Kendall C, Leonardi C, Hoffman PC, Combs DK. Intake and milk production of cows fed diets that differed in dietary neutral detergent fiber and neutral detergent fiber digestibility. J Dairy Sci. 2009;92(1):313–23.

    CAS  PubMed  Google Scholar 

  2. 2.

    Zhu W, Fu Y, Wang B, Wang C, Ye JA, Wu YM, Liu JX. Effects of dietary forage sources on rumen microbial protein synthesis and milk performance in early lactating dairy cows. J Dairy Sci. 2013;96(3):1727–34.

    CAS  PubMed  Google Scholar 

  3. 3.

    NRC: Nutrient Requirements of Dairy Cattle: Seventh Revised Edition. Washington. DC: The National Academies Press; 2001. p. 2001.

    Google Scholar 

  4. 4.

    Sun HZ, Wang DM, Wang B, Wang JK, Liu HY, Guan le L, Liu JX. Metabolomics of four biofluids from dairy cows: potential biomarkers for milk production and quality. J Proteome Res. 2015;14(2):1287–98.

    CAS  PubMed  Google Scholar 

  5. 5.

    Farmer H. Alfalfa imports: 122,500 tons from January to November and dropped 5.5% on year-on-year basis. [2019-12-31]. http://www.dairyfarmer.com.cn/nnyw_xjxm/2019-12-31/322763.chtml. (in Chinese).

  6. 6.

    Pang YZ, Liu YP, Li XJ, Wang KS, Yuan HR. Improving biodegradability and biogas production of corn Stover through sodium hydroxide solid state pretreatment. Energy Fuel. 2008;22(4):2761–6.

    CAS  Google Scholar 

  7. 7.

    Zhang R, Zhu W, Zhu W, Liu J, Mao S. Effect of dietary forage sources on rumen microbiota, rumen fermentation and biogenic amines in dairy cows. J Sci Food Agric. 2014;94(9):1886–95.

    CAS  PubMed  Google Scholar 

  8. 8.

    Wang B, Tu Y, Jiang LS, Liu JX. Effect of cereal straw and alfalfa hay diet on amino acid profile of gastrointestinal digesta in lactating dairy cows. J Anim Physiol Anim Nutr. 2018;102(2):421–8.

    CAS  Google Scholar 

  9. 9.

    Bergman EN, Brockman RP, Kaufman CF. Glucose metabolism in ruminants: comparison of whole-body turnover with production by gut, liver, and kidneys. Fed Proc. 1974;33(7):1849–54.

    CAS  PubMed  Google Scholar 

  10. 10.

    Aschenbach JR, Kristensen NB, Donkin SS, Hammon HM, Penner GB. Gluconeogenesis in dairy cows: the secret of making sweet milk from sour dough. IUBMB Life. 2010;62(12):869–77.

    CAS  PubMed  Google Scholar 

  11. 11.

    Loor JJ, Bionaz M, Drackley JK. Systems physiology in dairy cattle: nutritional genomics and beyond. Annu Rev Anim Biosci. 2013;1:365–92.

    PubMed  Google Scholar 

  12. 12.

    Bionaz M, Loor JJ. Ruminant metabolic systems biology: reconstruction and integration of transcriptome dynamics underlying functional responses of tissues to nutrition and physiological state. Gene Regul Syst Biol. 2012;6:109–25.

    Google Scholar 

  13. 13.

    Shahzad K, Bionaz M, Trevisi E, Bertoni G, Rodriguez-Zas SL, Loor JJ. Integrative analyses of hepatic differentially expressed genes and blood biomarkers during the peripartal period between dairy cows overfed or restricted-fed energy prepartum. PLoS One. 2014;9(6):e99757.

    PubMed  PubMed Central  Google Scholar 

  14. 14.

    Li L, Briskine R, Schaefer R, Schnable PS, Myers CL, Flagel LE, Springer NM, Muehlbauer GJ. Co-expression network analysis of duplicate genes in maize (Zea mays L.) reveals no subgenome bias. BMC Genomics. 2016;17(1):875.

    PubMed  PubMed Central  Google Scholar 

  15. 15.

    Huntington GB. Energy metabolism in the digestive tract and liver of cattle: influence of physiological state and nutrition. Reprod Nutr Dev. 1990;30(1):35–47.

    CAS  PubMed  Google Scholar 

  16. 16.

    Reynolds CK, Harmon DL, Cecava MJ. Absorption and delivery of nutrients for milk protein synthesis by portal-drained viscera. J Dairy Sci. 1994;77(9):2787–808.

    CAS  PubMed  Google Scholar 

  17. 17.

    Lomax MA, Baird GD. Blood flow and nutrient exchange across the liver and gut of the dairy cow. Effects of lactation and fasting. Br J Nutr. 1983;49(3):481–96.

    CAS  PubMed  Google Scholar 

  18. 18.

    Holter JB, West JW, Mcgilliard ML. Predicting ad libitum dry matter intake and yield of Holstein cows. J Dairy Sci. 1997;80(9):2188–99.

    CAS  PubMed  Google Scholar 

  19. 19.

    Krishnamoorthy L, Mahal LK. Glycomic analysis: an array of technologies. ACS Chem Biol. 2009;4(9):715–32.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Chen Y, Ding L, Ju H. In situ cellular glycan analysis. Acc Chem Res. 2018;51(4):890–9.

    CAS  PubMed  Google Scholar 

  21. 21.

    Griffin ME, Hsieh-Wilson LC. Glycan engineering for cell and developmental biology. Cell Chem Biol. 2016;23(1):108–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    D'Angelo G, Capasso S, Sticco L, Russo D. Glycosphingolipids: synthesis and functions. FEBS J. 2013;280(24):6338–53.

    CAS  PubMed  Google Scholar 

  23. 23.

    Wiegandt H: Glycolipids: Elsevier; 2011.

  24. 24.

    Zhang X, Kiechle FL. Review: Glycosphingolipids in health and disease. Ann Clin Lab Sci. 2004;34(1):3–13.

    PubMed  Google Scholar 

  25. 25.

    Uchimura K. Keratan sulfate: biosynthesis, structures, and biological functions. Methods Mol Biol (Clifton, NJ). 2015;1229:389–400.

    CAS  Google Scholar 

  26. 26.

    Pomin VH. Sulfated glycans in inflammation. Eur J Med Chem. 2015;92:353–69.

    CAS  PubMed  Google Scholar 

  27. 27.

    Vailati-Riboni M, Meier S, Burke CR, Kay JK, Mitchell MD, Walker CG, Crookenden MA, Heiser A, Rodriguez-Zas SL, Roche JR, et al. Prepartum body condition score and plane of nutrition affect the hepatic transcriptome during the transition period in grazing dairy cows. BMC Genomics. 2016;17(1):854.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Funderburgh JL. Keratan sulfate: structure, biosynthesis, and function. Glycobiology. 2000;10(10):951–8.

    CAS  PubMed  Google Scholar 

  29. 29.

    Davies Y, Lewis D, Fullwood NJ, Nieduszynski IA, Marcyniuk B, Albon J, Tullo A. Proteoglycans on normal and migrating human corneal endothelium. Exp Eye Res. 1999;68(3):303–11.

    CAS  PubMed  Google Scholar 

  30. 30.

    Caterson B, Melrose J. Keratan sulfate, a complex glycosaminoglycan with unique functional capability. Glycobiology. 2018;28(4):182–206.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Sun H, Wang B, Wang J, Liu H, Liu J. Biomarker and pathway analyses of urine metabolomics in dairy cows when corn Stover replaces alfalfa hay. J Anim Sci Biotechnol. 2016;7(1):49.

    PubMed  PubMed Central  Google Scholar 

  32. 32.

    Yoon J, Blumer A, Lee K. An algorithm for modularity analysis of directed and weighted biological networks based on edge-betweenness centrality. Bioinformatics (Oxford, England). 2006;22(24):3106–8.

    CAS  Google Scholar 

  33. 33.

    Contreras-Lopez O, Moyano TC, Soto DC, Gutierrez RA. Step-by-step construction of gene co-expression networks from high-throughput Arabidopsis RNA sequencing data. Methods Mol Biol. 2018;1761:275–301.

    CAS  PubMed  Google Scholar 

  34. 34.

    Tanaka KI, Xue Y, Nguyen-Yamamoto L, Morris JA, Kanazawa I, Sugimoto T, Wing SS, Richards JB, Goltzman D. FAM210A is a novel determinant of bone and muscle structure and strength. Proc Natl Acad Sci U S A. 2018;115(16):E3759–e3768.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Balazs A, Ruffert C, Hegyi E, Hritz I, Czako L, Takacs T, Szepes Z, Nemeth BC, Gervain J, Izbeki F, et al. Genetic analysis of the bicarbonate secreting anion exchanger SLC26A6 in chronic pancreatitis. Pancreatology. 2015;15(5):508–13.

    CAS  PubMed  Google Scholar 

  36. 36.

    Corbetta S, Eller-Vainicher C, Frigerio M, Valaperta R, Costa E, Vicentini L, Baccarelli A, Beck-Peccoz P, Spada A. Analysis of the 206M polymorphic variant of the SLC26A6 gene encoding a cl- oxalate transporter in patients with primary hyperparathyroidism. Eur J Endocrinol. 2009;160(2):283–8.

    CAS  PubMed  Google Scholar 

  37. 37.

    Monico CG, Weinstein A, Jiang Z, Rohlinger AL, Cogal AG, Bjornson BB, Olson JB, Bergstralh EJ, Milliner DS, Aronson PS. Phenotypic and functional analysis of human SLC26A6 variants in patients with familial hyperoxaluria and calcium oxalate nephrolithiasis. Am J Kidney Dis. 2008;52(6):1096–103.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Alper SL, Stewart AK, Chernova MN, Zolotarev AS, Clark JS, Vandorpe DH: Anion exchangers in flux: functional differences between human and mouse SLC26A6 polypeptides. Novartis Foundation Symp 2006, 273:107–119; discussion 119–125, 261–104.

  39. 39.

    Chernova MN, Jiang L, Friedman DJ, Darman RB, Lohi H, Kere J, Vandorpe DH, Alper SL. Functional comparison of mouse slc26a6 anion exchanger with human SLC26A6 polypeptide variants: differences in anion selectivity, regulation, and electrogenicity. J Biol Chem. 2005;280(9):8564–80.

    CAS  PubMed  Google Scholar 

  40. 40.

    Ha JY, Kim JS, Kang YH, Bok E, Kim YS, Son JH. Tnfaip8 l1/Oxi-beta binds to FBXW5, increasing autophagy through activation of TSC2 in a Parkinson's disease model. J Neurochem. 2014;129(3):527–38.

    CAS  PubMed  Google Scholar 

  41. 41.

    Hu J, Zacharek S, He YJ, Lee H, Shumway S, Duronio RJ, Xiong Y. WD40 protein FBW5 promotes ubiquitination of tumor suppressor TSC2 by DDB1-CUL4-ROC1 ligase. Genes Dev. 2008;22(7):866–71.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Golob-Schwarzl N, Wodlej C, Kleinegger F, Gogg-Kamerer M, Birkl-Toeglhofer AM, Petzold J, Aigelsreiter A, Thalhammer M, Park YN, Haybaeck J. Eukaryotic translation initiation factor 6 overexpression plays a major role in the translational control of gallbladder cancer. J Cancer Res Clin Oncol. 2019;145(11):2699–711.

    CAS  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Pinzaglia M, Montaldo C, Polinari D, Simone M, La Teana A, Tripodi M, Mancone C, Londei P, Benelli D. EIF6 over-expression increases the motility and invasiveness of cancer cells by modulating the expression of a critical subset of membrane-bound proteins. BMC Cancer. 2015;15:131.

    PubMed  PubMed Central  Google Scholar 

  44. 44.

    Kim SE, Hinoue T, Kim MS, Sohn KJ, Cho RC, Weisenberger DJ, Laird PW, Kim YI. Effects of folylpolyglutamate synthase modulation on global and gene-specific DNA methylation and gene expression in human colon and breast cancer cells. J Nutr Biochem. 2016;29:27–35.

    CAS  PubMed  Google Scholar 

  45. 45.

    Freemantle SJ, Taylor SM, Krystal G, Moran RG. Upstream organization of and multiple transcripts from the human folylpoly-gamma-glutamate synthetase gene. J Biol Chem. 1995;270(16):9579–84.

    CAS  PubMed  Google Scholar 

  46. 46.

    Kim YI. Folate and colorectal cancer: an evidence-based critical review. Mol Nutr Food Res. 2007;51(3):267–92.

    CAS  PubMed  Google Scholar 

  47. 47.

    Shane BJFih, disease: Folate chemistry and metabolism. 1995:1–22.

  48. 48.

    Zhou XQ, Zhang YD, Zhao M, Zhang T, Zhu D, Bu DP, Wang JQ. Effect of dietary energy source and level on nutrient digestibility, rumen microbial protein synthesis, and milk performance in lactating dairy cows. J Dairy Sci. 2015;98(10):7209–17.

    CAS  PubMed  Google Scholar 

  49. 49.

    Gao ST, Ma L, Zhou Z, Zhou ZK, Baumgard LH, Jiang D, Bionaz M, Bu DP. Heat stress negatively affects the transcriptome related to overall metabolism and milk protein synthesis in mammary tissue of midlactating dairy cows. Physiol Genomics. 2019;51(8):400–9.

    PubMed  Google Scholar 

  50. 50.

    Bu D, Bionaz M, Wang M, Nan X, Ma L, Wang J. Transcriptome difference and potential crosstalk between liver and mammary tissue in mid-lactation primiparous dairy cows. PLoS One. 2017;12(3):e0173082.

    PubMed  PubMed Central  Google Scholar 

  51. 51.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Love MI. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.

    Google Scholar 

  54. 54.

    Bionaz M, Periasamy K, Rodriguez-Zas SL, Hurley WL, Loor JJ. A novel dynamic impact approach (DIA) for functional analysis of time-course omics studies: validation using the bovine mammary transcriptome. PLoS One. 2012;7(3):e32455.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    PubMed  Google Scholar 

  56. 56.

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

    CAS  Google Scholar 

Download references

Acknowledgements

We give especial thanks to our former and current students who have contributed the projects.

Funding

The design and sample collection and detection of this study were supported by the National Key Research and Development Program of China (2018YFD0501600). The Scientific Research Project for Major Achievements of The Agricultural Science and Technology Innovation Program (ASTIP) (CAAS-ZDXT2019004, ASTIP-IAS07–1, CAAS-XTCX2016011–01) played critical roles in data analysis and interpretation of this study. Manuscript writing and revision of this study were supported by Beijing Dairy Industry Innovation Team (BAIC06–2020).

Author information

Affiliations

Authors

Contributions

DB and JW designed the study; YZ and LM conducted the animal trial and collected samples; DB and LM prepared libraries; SG wrote the manuscript and analyzed data. JL revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to D. P. Bu.

Ethics declarations

Ethics approval and consent to participate

All animal care and procedures were approved by the Animal Welfare and Ethical Committee of Institute of Animal Science, Chinese Academy of Agricultural Sciences (No. IAS20180115). The use of animals in the present study was consented in written form by Beijing ZhongDi Dairy Farm. The whole experiment was conducted in strict accordance with the Directions for Caring of Experimental Animals from the Institute of Animal Science, Chinese Academy of Agricultural Sciences.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1.

Dataset of differently expressed genes (DEG) of CS vs. MF.

Additional file 2.

Output of Dynamic Impact Approach (DIA) analysis of the DEG.

Additional file 3.

Output of Database for Annotation, Visualization and Integrated Discovery (DAVID) analysis of the DEG.

Additional file 4.

The table containing the statistically significant correlations across the whole data set for every pair of DEGs.

Additional file 5.

Co-expression network statistics including correlation of the DEGs, the degree and betweenness centrality of the nodes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Gao, S.T., Ma, L., Zhang, Y.D. et al. Hepatic transcriptome perturbations in dairy cows fed different forage resources. BMC Genomics 22, 35 (2021). https://doi.org/10.1186/s12864-020-07332-0

Download citation

Keywords

  • Liver transcriptome
  • Forage resources
  • RNAseq
  • Corn Stover