Hepatic transcriptome perturbations in dairy cows fed different forage resources
BMC Genomics volume 22, Article number: 35 (2021)
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.
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”.
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.
Forage is the largest component of the diet of lactating cows and could affect dry matter intake (DMI)  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 . Increasing RUP by 1% can improve milk production by 1 kg . 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 .
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 . In 2019, more than 1.2 million tons of alfalfa hay were imported to cover the shortage of high-quality forage in China . At the same time, it is estimated that more than 10 million tons of corn stover are generated annually in China . 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 .
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 . 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 . 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 .
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 . 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.
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.
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).
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).
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).
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 . 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.
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).
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 . Metabolic function and, thus, energy metabolism of liver responds to a variety of environmental stimuli including fasting or level of feed intake , diet composition and productive (physiological) state . 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 , 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 . 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) , 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 . 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 , and mediate a wide variety of biological processes including cell growth and differentiation, cell−cell communication, immune response, pathogen interaction, and intracellular signaling events . 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) . 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 , and are classified based on their carbohydrate structure into six major series in vertebrates including gangliosides, lacto-, neolacto-, muco-, isoglobo-, and globo-series GSL . 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 . 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 . 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 . 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 . 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 . 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 . Thus, it was suggested that overfeeding of fatter cows may decrease the synthesis of anti-inflammatory compounds and consequently induce some detrimental effects . 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 . 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 . 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 . These metabolites were involved in glycine, serine, and threonine metabolism; tyrosine metabolism; and phenylalanine metabolism . 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 . 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 . 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 . 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 . Thus, both degree and betweenness centrality are measures of the function of a node in network connectivity . 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 . 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 . 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 . Ha et al. (2014) demonstrated that intracellular accumulation of TSC2 inhibits the activity of mTOR and increase autophagy . 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 , while EIF6 over-expression increases the motility and invasiveness of cancer cells . 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.
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.
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) . 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.
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 . 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 . 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 .
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 . 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 . 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 . 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) .
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” (220.127.116.11) was used to calculate network statistics (degree and betweenness centrality for each node) as described by Contreras-Lopez et al. (2018) . 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.
Acid detergent fiber
Armadillo repeat containing X-linked 2
Database for Annotation, Visualization and Integrated Discovery
Damage specific DNA binding protein 1
Differently expressed genes
Dynamic Impact Approach
Dry matter intake
Eukaryotic translation initiation factor 6
Family with sequence similarity 210 member A
F-box and WD repeat domain containing 5
False discovery rates
Gene Ontology biological process
Gene Ontology cellular components
Gene Ontology molecular function
Kyoto Encyclopedia of Genes and Genomes
Human gallbladder cancer
Neutral detergent fiber
Rumen degradable protein
Rotamase CYP 1
Rumen undegradable proteins
Solute carrier family 26 member 6
Human extrahepatic bile duct carcinoma cell line
TSC complex subunit 2
Zinc finger and SCAN domain containing 10
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.
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.
NRC: Nutrient Requirements of Dairy Cattle: Seventh Revised Edition. Washington. DC: The National Academies Press; 2001. p. 2001.
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.
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).
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.
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.
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.
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.
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.
Loor JJ, Bionaz M, Drackley JK. Systems physiology in dairy cattle: nutritional genomics and beyond. Annu Rev Anim Biosci. 2013;1:365–92.
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.
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.
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.
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.
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.
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.
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.
Krishnamoorthy L, Mahal LK. Glycomic analysis: an array of technologies. ACS Chem Biol. 2009;4(9):715–32.
Chen Y, Ding L, Ju H. In situ cellular glycan analysis. Acc Chem Res. 2018;51(4):890–9.
Griffin ME, Hsieh-Wilson LC. Glycan engineering for cell and developmental biology. Cell Chem Biol. 2016;23(1):108–21.
D'Angelo G, Capasso S, Sticco L, Russo D. Glycosphingolipids: synthesis and functions. FEBS J. 2013;280(24):6338–53.
Wiegandt H: Glycolipids: Elsevier; 2011.
Zhang X, Kiechle FL. Review: Glycosphingolipids in health and disease. Ann Clin Lab Sci. 2004;34(1):3–13.
Uchimura K. Keratan sulfate: biosynthesis, structures, and biological functions. Methods Mol Biol (Clifton, NJ). 2015;1229:389–400.
Pomin VH. Sulfated glycans in inflammation. Eur J Med Chem. 2015;92:353–69.
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.
Funderburgh JL. Keratan sulfate: structure, biosynthesis, and function. Glycobiology. 2000;10(10):951–8.
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.
Caterson B, Melrose J. Keratan sulfate, a complex glycosaminoglycan with unique functional capability. Glycobiology. 2018;28(4):182–206.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Kim YI. Folate and colorectal cancer: an evidence-based critical review. Mol Nutr Food Res. 2007;51(3):267–92.
Shane BJFih, disease: Folate chemistry and metabolism. 1995:1–22.
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.
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.
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.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.
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.
Love MI. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.
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.
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.
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.
We give especial thanks to our former and current students who have contributed the projects.
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).
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Dataset of differently expressed genes (DEG) of CS vs. MF.
Output of Dynamic Impact Approach (DIA) analysis of the DEG.
Output of Database for Annotation, Visualization and Integrated Discovery (DAVID) analysis of the DEG.
The table containing the statistically significant correlations across the whole data set for every pair of DEGs.
Co-expression network statistics including correlation of the DEGs, the degree and betweenness centrality of the nodes.
About this article
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
- Liver transcriptome
- Forage resources
- Corn Stover