Assembly and Analysis of Changes in Transcriptomes of Dairy Cattle Rumen Epithelia during Lactation and Dry Periods ()
1. Introduction
Lactation in dairy cattle is inextricably coupled with an increased nutrient requirement to support milk synthesis. Therefore, along with massive increases in feed intake the underlying metabolism of the cow must adapt to lactation-associated challenges and requires major functional adjustments of the rumen and the entire digestive system. Dairy cattle are typically fed a ration with low to moderate fermentability during the early dry period (non-lactating) and with increased fermentability as the cows near parturition. After parturition of course, diet fermentability is further increased to meet the energy demands arising from the onset of lactation [1] . Therefore, the rumen epithelium is faced with an ever-changing ruminal environment and thus, metabolic reprogramming of the rumen epithelial cells is likely characterized by alterations in the regulation of pathways in order to provide the cow with vital energy and metabolites to support the physiologic state.
As the consequence of feeding highly fermentable diets to ruminants to increase energy intake during lactation period, an increase in short-chain fatty acid (SCFAs; acetate, propionate, and butyrate) production and a reduced ruminal pH create challenges to metabolism and the regulation of intracellular pH homeostasis of the ruminal epithelia [2] . The epithelia respond to these challenges in a coordinated manner. An enlargement of the absorptive surface area is well documented [2] , but emerging evidence indicates that changes in epithelial cell functions may be the initial response [3] . The SCFAs are important nutrients in ruminants. SCFAs are produced during the microbial fermentation of dietary fiber in the gastrointestinal tract and are directly absorbed at the site of production and oxidized for cell energy production and use [4] . The rumen is the site of microbial fermentation of dietary fiber to SCFA, and it functions as the focal point of energy digestion and absorption of the ruminant gastrointestinal tract. SCFAs are absorbed by the rumen epithelium and account for over 70% of the total energy requirements available for the cow [3] [5] [6] . It is very reasonable to assume that rumen epithelial adaptation to the challenges of metabolism and the high nutrient requirement during lactation is a primary trigger for alterations in genetic activities, and gene expression in particular. However, little was known regarding the rumen epithelial transcriptome, changes induced by lactation, as well as functional adaptation of the rumen epithelium after the transition until recently. A published report using a micro array to characterize the adaptation of the rumen epithelium during transition, mRNA profiling identified a total of 1011 (−3 vs. +1 wk) and 729 (−3 vs. +6 wk) differentially expressed genes [7] . In another report, a total of 35 genes associated with metabolism, transport, inflammation, and signaling were evaluated by quantitative reverse transcription-PCR and to investigate the changes in rumen microbiota, gene expression of the ruminal epithelium, and blood biomarkers of metabolism and inflammation during the transition period. The results indicated that alterations in ruminal epithelium gene expression could be driven by nutrient intake-induced changes in microbes; microbial metabolism; and the systemic metabolic, hormonal, and immune changes [8] .
High-throughput sequencing of RNA (RNA-Seq) is an efficient way of mapping and quantifying transcriptome and was developed to help analyze global gene expression in cells and tissues. Unlike microarray or PCR based transcriptomic analysis, RNA-Seq does not use a pre-defined set of probes to capture and quantify RNA sequences within a sample. Therefore RNA-Seq can be used to investigate both known and novel transcripts. A further advantage of RNA-Seq is a large dynamic range of expression levels. In this report, we used next-generation sequencing technology to assemble and to profile and compare the transcriptome of cattle rumen epithelial tissues from cattle in either their dry or lactation periods. Our data indicated substantial dynamics of gene expression during the dry and lactation periods, thereby supporting the importance of assessment of rumen epithelial metabolism by molecular methods. This study also provides crucial information to support further exploration of the functional adaptation of epithelial cells to the challenge of the nutritional requirements for lactation.
2. Materials and Methods
2.1. Animals, Experimental Treatments
All animal care and handling were conducted according to the guidelines approved by the USDA Beltsville Area Institutional Animal Care Committee. Four first lactation cows (about 3 years old) were selected randomly from the Beltsville Dairy Herd. Cows were ruminally cannulated as previously described (Li et al., 2012, Baldwin et al., 2012) during second week of the dry period. Samples of ruminal tissue were collected from cows during their second lactation (approximate 200 days in milk) and Dry cow samples collected during the dry period between their 2nd and 3rd lactation. Mid-lactation dairy cows were fed ad libitum the standard lactation rations with approximately 50% corn silage and 50% concentrate on a dry matter basis (53.3% DM, 16.4% CP, 65.5% Degradable protein and an NEL of 1.7 Mcal/kg). Similarly, Dry cows were fed a dry cow ration comprised of 17% concentrate with 29% wheat straw, 22.7% corn silage, and 31% haylage on a dry matter basis (52.6% DM, 13.2% CP, 70% Degradable protein with NEL 1.5 Mcal/kg). Rumen epithelial samples from the dorsal sac of the rumen immediately adjacent to the reticulo-rumen orifice were serially collected via grab biopsy by pinching with a gloved hand through rumen fistulae.
2.2. RNA Isolation and Sequencing
Methods for RNA isolation and sequencing were previously described in our earlier publication [9] . Briefly, total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) followed by DNase digestion and Qiagen RNeasy column purification (Qiagen, Valencia, CA, USA). The RNA integrity was verified using an Agilent Bioanalyzer 2100 (Agilent, Palo Alto, CA, USA). High-quality RNA (RNA Integrity number or RIN > 9.0) was processed using an Illumina TruSeq RNA sample prep kit following the manufacturer’s instruction (Illumina, San Diego, CA, USA). After quality control procedures, individual RNA-seq libraries were then pooled based on their respective 6-bp adaptors and sequenced at a 50bp/sequence and read using an Illumina HiSeq 2000 sequencer platform, as described previously [10] . Approximately 55 million reads per sample (mean ± sd = 55.09 million ± 10.42 million) were generated.
2.3. Data Analysis and Bioinformatics
The CLC Genomics Workbench (v10, Qiagen) was used for RNA-sequencing data analysis. Raw sequence reads were first checked using a quality control pipeline. Nucleotides of each raw read were scanned for low quality and trimmed. Trimmed reads were aligned to the bovine reference genome (BosTau_UMD3.1). Gene expression levels were normalized as reads per kilobase of exon model per million mapped reads (RPKM) using the CLC transcriptomic analysis tool. The top 200 highly expressed genes selected by the sum of RPKM values in all libraries were used for heat map analysis with hierarchical clustering module in the CLC Genomics Workbench (Setting: sample and feature clustering, normalized expression value, Euclidean distance and complete linkage).
Differentially expressed genes in the transcriptome were further analyzed using a gene ontology (GO) enrichment analysis module in the CLC Genomics Workbench. Enrichment of certain GO terms was determined based on Fisher’s exact test. A multiple correction control (permutation to control false discovery rate) was implemented to set up the threshold to obtain the lists of significantly over-represented GO terms.
The molecular processes, molecular functions, and genetic networks were further evaluated by analyzing differentially expressed genes using Ingenuity Pathways Analysis (IPA, Ingenuity® Systems, and http://www.ingenuity.com). IPA is a software application that enables biologists to identify the biological mechanisms, pathways, and functions most relevant to their experimental datasets or genes of interest [11] [12] [13] [14] [15] .
2.4. Canonical Pathway Analysis of Data Sets
Analysis of canonical pathways identified the pathways from the IPA library of canonical pathways that were most significant to the data set. Genes from the data set that were associated with a canonical pathway in the Ingenuity Pathways Knowledge Base were considered for the analysis. The significance of the association between the data set and the canonical pathway was measured in two ways: one was to use a ratio of the number of genes from the data set that map to the pathway divided by the total number of genes that map to the canonical pathway was displayed. The other was to use Fischer’s exact test to calculate a p-value that determines the probability that the association between the genes in the dataset and the canonical pathway was explained by chance alone.
2.5. Pathway Analysis and Network Generation
The dataset containing gene identifiers and corresponding expression values was uploaded into the Ingenuity Pathways Analysis application. Each gene identifier was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base. These genes, called focus genes, were overlaid onto a global molecular network developed from information contained in the Ingenuity Pathways Knowledge Base. Networks of these focus genes were then algorithmically generated based on their connectivity.
3. Results
3.1. RNA Sequencing Data Statistics and the Rumen Epithelial Transcriptome Profiling
The molecular mechanisms underlying differential metabolism between lactating and dry cows were elucidated by creating eight cDNA libraries from four independent biological samples (rumen papillae) of lactating and dry cows (n = 4 per group) and sequenced using the Illumina HiSeq 2000ä sequencing platform. About 55 million (55.09 million ± 10.42 million, mean ± SD) high-quality reads were retrieved for each sample. Cleaned reads were de novo mapped to the reference genome (BosTau_UMD3.1). In total, 96.96% ± 1.35% (mean ± SD) of the sequence reads were mapped to the reference genome. The details of the mapping results are listed in Table 1. The mapping results indicate a high quality sequencing and transcriptomic assembly. A total 24,956 genes were detected at least once in one of the rumen epithelial samples. The number of genes expressed per sample varied from 15,184 to 16,687 (15765 ± 610, mean ± SD).
![]()
Table 1. RNA sequencing mapping results.
3.2. Global Gene Expression Dynamics between Dry and Lactation Periods
The PCA and volcano plot of differentially expressed genes (Figure 1(a) and Figure 1(b)) revealed significant dynamic changes of global gene expression between rumen epithelia from dry cows and lactating cattle. The PCA plot shows that gene expression profiles of rumen epithelia from dry and lactation cattle fall into two very distinct clusters (Figure 1(a)). Unsupervised hierarchical clustering of top 100 differentially expressed genes (Figure 2) also indicates the distinctive gene expression patterns between lactating and dry cattle. We found 5036 genes to be significantly differently expressed (DE) with FDR p-value < 0.05 and
![]()
Figure 1. (a) The PCA plot shows that gene expression profiles of rumen epithelia from dry and lactation cattle fall into two very distinct clusters. (b) Volcano plots: the differentially expressed genes during dry and lactation periods.
![]()
Figure 2. Heat map. Hierarchical clustering of selected 100 genes significantly differentially expressed genes during day and lactation periods in rumen epithelium. The expression scale is log2 based.
at least 5-fold difference (FPKM). In lactating cattle, 2107 genes were up-regulated and 2929 genes were down-regulated in comparison to rumen epithelia from dry cattle. All significant differentially expressed genes are listed in Supplementary Table 1. The top 10 up-regulated unigenes and top 10 down-regulated unigenes in rumen epithelial tissue in lactation cattle are listed in Table 2.
3.3. Gene Ontology (GO) Enrichment Analysis of Genes Impacted during Lactation in Cattle
Gene ontology (GO) enrichment analysis was performed for further exploration of the potential functions of the differentially expressed genes in regulating cattle rumen function during the lactation period. The GO enrichment of significant differentially expressed genes was categorized into 24 groups in the biological process (FDR p-value < 0.05, Table 2). The most significantly enriched GO term in the biological process was a carbohydrate derivative metabolic process and was followed by nucleoside metabolic processes (Table 3). Interestingly, most
![]()
Table 2. The top 10 up-regulated and top 10 down regulated unigenes in rumen epithelial tissue of lactation cattle.
![]()
Table 3. Gene Ontology (GO) enrichment analysis in biological process shows GO terms significant changed in lactation cow.
GO terms in the 24 groups were related to various metabolic processes. Additionally, the GO term in molecular function analysis showed only two groups: catalytic activity and transferase activity with FDR p-value < 0.05 (Table 4). These results indicated that the transcriptomic changes in rumen epithelia of lactating cattle were triggered by rumen epithelial adaptation to the challenges of metabolism and high nutrient requirements during the lactation. The GO term in the cellular component analysis was also determined (Table 5).
3.4. Function and Regulatory Gene Network Analysis of Differentially Expressed Genes Using IPA
We looked further into the biological function and differentially expressed genes in lactating cattle by performing function and pathway analysis using IPA. Analysis of differentially expressed genes in lactating cattle revealed the most significantly perturbed biological functions in the dataset. The differentially expressed genes were first analyzed as a whole dataset. The top 14 molecular and cellular functions, as determined by p-value, are shown in Figure 3(a). These functions, such as cell morphology, cellular function, and maintenance, characterized the mechanisms related to the adaptation of rumen epithelia to the physiological changes during lactation, which were represented in the physiological function changes revealed by IPA analysis (Figure 3(b)). We also specifically analyzed the up-regulated genes and down-regulated genes separately in the dataset (Figure 4(a) and Figure 4(b)). Those functions could exemplify the most significantly enhanced or inhibited biological functions in the rumen epithelia of lactating cattle.
![]()
Table 4. Gene Ontology (GO) enrichment analysis in molecular function shows GO terms significant changed in lactation cow.
![]()
Table 5. Gene Ontology (GO) enrichment analysis in cellular component.
![]()
(a) ![]()
(b)
Figure 3. Global functional analyses. (a) Molecular and cellular function analysis of the differentially expressed genes; (b) Physiological functions of the differentially expressed genes. The significance value associated with a function in global analysis is a measure for how likely it is that genes from the dataset file under investigation participate in that function. The significance is expressed as a p value that is calculated using the right-tailed Fisher’s exact test.
3.5. IPA Analysis Unveiled Potential Upstream Regulators in the Dataset
Potential upstream regulators may have an involvement in the activities of the perturbed genes during lactation. Those upstream regulators can be divided into different groups depending on their molecular types. The two most relevant groups are transcription regulators and ligand-dependent nuclear receptors. IPA upstream regulator analysis identified 42 transcription regulators in the dataset (Listed in Supplementary Table 2). Most of these transcription regulators in the
![]()
(a)![]()
(b)
Figure 4. Global functional analyses. (a) Molecular and cellular functions of up-regulated genes; (b) Molecular and cellular functions of down regulated genes.
dataset were predicted in activated states, while seven of the transcription regulators (COMMD3-BMI1, KDM5A, ASXL1, FBXW7, SOX3, RBPJ, HAND1 and MED13) were in the inhibited state. Another major category of upstream regulators was classified as ligand-dependent nuclear receptors. Remarkably, three members of the Peroxisome Proliferator Activated Receptor (PPAR) family, PPARA (PPAR-α), PPARD (PPAR-δ) and PPARG (PPAR-γ), were identified as the activated upstream regulators in this group (Table 6).
![]()
Table 6. Upstream regulators-ligand-dependent nuclear receptors.
3.6. The Biologically Relevant Regulatory Gene Networks
The IPA also exposed the biologically relevant regulatory networks in the up-regulated gene dataset. The full list of networks is presented in Supplementary Table 3. The top network in the up-regulated gene group (Figure 5) indicates that the function of molecule transport is increased in the rumen epithelia of lactating cattle. We identified additional activated molecular transporters during lactation by searching for the dataset of up-regulated genes, and we discovered more than 60 molecule transporter genes. Those transporter genes and their functional annotations are listed in Table 7.
Another functionally relevant gene network is the network centered on the PPARA gene (Figure 6). The major functions of this network are lipid metabolism, molecule transport, and small molecule biochemistry. The PPARA gene was identified as an upstream regulator of gene expression (Table 6) and was in the center of this network.
4. Discussion
Many of the recent advances in biology have been driven by genome sequence information. However, the next significant challenge is to understand the complexity of genome information to predict the phenotypes, to annotate functional elements encoded in the genome, and to understand the fundamental consequences of the genome, which are dictated by transcriptomics. In this report, next generation sequencing and transcriptomic profiling were used to investigate
![]()
Figure 5. The top network in the up-regulated gene group indicates that the function of molecule transport is increased in the rumen epithelia of lactating cattle. The note color indicates the expression level of the genes. Notes and edges are displayed with various shapes and labels that present the functional class of genes and the nature of the relationship between the notes, respectively. The red notes indicate up-regulated genes.
![]()
Figure 6. A functionally relevant gene network centered on the PPARA gene indicates that PPAR family proteins play various roles in the functional adaptation of rumen metabolism to meet the nutritional demands of lactation.
![]()
Table 7. Up-regulated molecular transporters.
the molecular basis of the adaptation of rumen epithelia of dairy cattle to provide the increased nutrient requirements for milk synthesis. We assembled the transcriptome and compared gene expression patterns in rumen epithelial tissue of dairy cattle in both dry and lactation periods. PCA analysis and volcano plots revealed that gene expression profiles of the rumen epithelia from dry and lactating cattle fell into two very distinct clusters.
The very distinctive rumen epithelial transcriptome reflects the morphological and functional differences of rumens from the dairy cattle in the dry period versus the lactation period. Notably, the four top up-regulated unigenes were MYH1 and MYH7 (myosin heavy chain 1 and 7) and related genes MYBPC1 and MYL2. Myosin is a major contractile protein that converts chemical energy into mechanical energy through the hydrolysis of ATP. A multigene family encodes myosin heavy chains and they are expressed in different stages of animal development and in different fiber types. The MYH gene family shows expression that is spatially and temporally regulated during development [16] . Up-regulated expression of these genes may indicate an involvement in rumen morphological changes during lactation and increased cellular movement of rumen papilla [7] .
The rumen epithelium is metabolically active and is the greatest consumer of energy of the total viscera [2] [17] . Information generated from this work supports our understanding of the functional adaptation of rumen metabolism to meet the nutritional demands of lactation. Milk production in dairy cattle substantially increases their dietary nutrient and energy requirements to support milk synthesis. Therefore, the metabolism of dairy cattle must adapt to milk production-associated challenges that require major functional adjustments of the rumen and digestive system. Increased nutrient requirements for milk synthesis trigger greater dry matter intake and lead to morphological changes in rumen of dairy cattle, such as increases in rumen papilla surface area and result in a net increase in absorption rate of short-chain fatty acids [18] [19] [20] after calving. Our data indicated that rumen morphological structure and expression of genes in dairy cattle during lactation are correlated with metabolic adaptation to the demanding nutrient requirements [21] . The effects of enhanced dietary plane of nutrition on the gene expression pattern of ruminal epithelial tissue of young Holstein calves also indicated that enhanced dietary nutrition could elicit a strong transcriptomic response in the ruminal epithelial tissue [22] .
Gene ontology (GO) enrichment analysis in biological processes shows GO terms significant changed in dairy cattle during lactation. Most of the genes perturbed in the rumen epithelial tissues of dairy cattle during lactation are related to the various metabolic processes in GO terms of biological processes. Functional analysis performed by IPA also yielded results consistent with the finding that these functions, such as cell morphology, cellular function, and maintenance, characterize the mechanisms related to the adaptation of rumen epithelia to the physiological changes during the lactation.
Members of the PPAR (Peroxisome Proliferator Activated Receptor) family are known to bind to metabolites such as fatty acids and affect the transcription of many genes [21] . PPARα is the primary transcription factor responsible for the induction of the majority of genes necessary for fatty acid transport, uptake, and oxidation, as well as ketone body biosynthesis and import [23] [24] . PPARα, a multi-functional protein and one of the key transcription factors in the metabolic network, participates in signaling driven by the main nutrient sensors, such as AMP-activated protein kinase (AMPK), PPARγ coactivator1α (PGC-1α) [25] . PPARα is also a master regulator for tissue development and is responsive to weaning in rumen tissue [26] . The transcriptional activities of PPAR family are dynamic, yet highly coordinated. PPARα may synchronize the dynamic regulation of related gene expression [2] . An additional level of complexity in the PPARA regulatory mechanism arises from its ability to interact with various partners, which ultimately determines the metabolic phenotype [27] . Our study indicates that PPAR family proteins play various roles in the functional adaptation of rumen metabolism to meet the nutritional demands of lactation. Therefore, exploring the regulatory roles of the PPAR family of transcription factors and the functions of their target gene groups may be an important step in molecular research into the metabolic functions of rumen epithelia that provide the increased nutrient requirements for milk synthesis.
In addition to genes related to cellular movement and tissues morphology, those involved in molecular transport have been identified as differentially expressed genes during adaptation to lactation. Extreme changes in the rumen environment occur during lactation due to a highly fermentable diet. The increased supply of SCFAs is among the most important changes [2] and as discussed earlier the ruminal epithelium has an enormous capacity for the absorption of SCFAs. This not only delivers metabolic energy to the animal but it is also an essential regulatory mechanism that stabilizes the intraruminal milieu [3] . However, the influx of SCFAs may cause severe cytosolic acidosis. The ruminal epithelium can buffer protons by various mechanisms, such as via a Na+/H+ exchanger, a bicarbonate importing system, or an H+/monocarboxylate cotransporter (MCT) [3] . In the dataset and the functional network (Figure 5), ATP1B1 (ATPase Na+/K+ transporting subunit beta 1) plays an important role in this network. In addition, NR3C1, NR3C2 (nuclear receptor subfamily 3 group C member 1 and C member 2), and ACSL4 are also major components of this network.
The protein encoded by the ACSL4 gene is an isozyme of the long-chain fatty-acid-coenzyme A ligase family. All isozymes of this family transform free long-chain fatty acids into fatty acyl-CoA esters, and thus play a key role in lipid biosynthesis and fatty acid degradation. This protein also functions as a molecular transporter [28] . Increased transport activity in the rumen epithelial tissue in lactating cattle represents an epithelial adaptation associated with the substantial diet change to meet the nutrient demand of lactation. Increased transport activities also reflect the increased availability of energy sources, such as SCFAs. Increases in the absorption rates of SCFAs are well documented when animals are adapting to energy-rich diets [2] [3] .
The nutritional strategies for maximizing milk production and health in lactation are continually developing [7] . Traditionally, dairy cattle were fed increasingly more energy rich diets during the pre-parturition period to adapt the rumen, in order to facilitate an easier transition to the higher energy density of diets fed to lactating cows. The data presented in this report from implementation of high-throughput transcriptomic data analysis are coordinated well with rumen function. Transcriptomics profiling and comparison reveal that the adaption to lactation was marked by extensive changes in gene expression related to metabolism in rumen epithelial tissue. Up-stream regulators, such as PPARA, and genes of molecular transporters are the points of focus revealed in this report. Our report not only confirmed earlier studies but also provides additional information using robust technology and bioinformatics analysis on the distinct features of the rumen at two important physiologic stages of dairy cattle. Information generated from the work enhances our understanding of the functional adaptation of rumen metabolism to meet the nutritional demands of lactation and promotes future advances in ruminant nutrition and feeding.
Acknowledgements
Mention of a product, reagent or source does not constitute an endorsement by USDA to the exclusion of other products or services that perform a comparable function. The US Department of Agriculture is an equal opportunity provider and employer.
Author Contributions
RLB, and CJL conceived and performed animal experiment, RWL prepared RNA library and RNA sequencing. CJL performed RNA-seq data analysis and wrote the manuscript. RLB, RWL and CJL all read and proofed the manuscript.
Competing Financial Interests
Authors declare no conflict of interest.