Transcriptome Analysis of Bovine Rumen Tissue in Three Developmental Stages.

Rumen development is a crucial physiological challenge for ruminants. However, the molecular mechanism regulating rumen development has not been clearly elucidated. In this study, we investigated genes involved in rumen development in 13 rumen tissues from three developmental stages (birth, youth, and adult) using RNA sequencing. We identified that 6,048 genes were differentially expressed among three developmental stages. Using weighted correlation network analysis, we found that 12 modules were significantly associated with developmental stages. Functional annotation and protein-protein interaction (PPI) network analysis revealed that CCNB1, CCNB2, IGF1, IGF2, HMGCL, BDH1, ACAT1, HMGCS2, and CREBBP involved in rumen development. Integrated transcriptome with GWAS information of carcass weight (CW), stomach weight (SW), marbling score (MS), backfat thickness (BFT), ribeye area (REA), and lean meat weight (LMW), we found that upregulated DEGs (fold change 0∼1) in birth-youth comparison were significantly enriched with GWAS signals of MS, downregulated DEGs (fold change >3) were significantly enriched with GWAS signals of SW, and fold change 0∼1 up/downregulated DEGs in birth-adult comparison were significantly enriched with GWAS signals of CW, LMW, REA, and BFT. Furthermore, we found that GWAS signals for CW, LMW, and REA were enriched in turquoise module, and GWAS signals for CW was enriched in lightgreen module. Our study provides novel insights into the molecular mechanism underlying rumen development in cattle and highlights an integrative analysis for illustrating the genetic architecture of beef complex traits.

Introduction

Rumen development is an important physiological challenge for young ruminants (Lin et al., 2019). At birth, rumen is incompletely developed without high ketogenic capacity (Warner et al., 1956; Lin et al., 2020). With the rumen development and microbial colonization, calves began to transit from a milk-based diet to a grain and forage-based diet (Jami et al., 2014; Shabat et al., 2016). The development of the rumen during this transition involves three simultaneous processes ( ). First is the physical development of rumen including growth in rumen volume and papilla (Reynolds et al., 2004). Second, microbial communities were established and colonized (Fouts et al., 2012; Rey et al., 2014), which is highly correlated with important carcass traits, such as marbling score (MS), adjusted 12th rib fat thickness, longissimus lipid content, and carcass yield grade ( ). Meanwhile, there is a functional development of fermentation capacity and enzyme activity in the rumen lumen and epimural layers (Rey et al., 2012). This whole process is not instantaneous, and the ability to ferment solid diet acquired at least 1 week after weaning (Quigley et al., 1985). However, little is known about transcriptome characteristics during these morphological changes.

Previous studies focused on the effects of nutrient, diet composition, or feeding strategy on rumen development. Encouraging calves to consume dry feedstuffs at an early age will accelerate rumen development, allowing more efficient body growth and development at maturity (Coverdale et al., 2004; Norouzian et al., 2011). Commonly rumen microorganisms proliferate and produce energy in the form of volatile fatty acids (VFAs), primarily propionate and butyrate (Reddy et al., 2017). VFAs, particularly butyrate, can stimulate papilla growth, accelerate rumen motility, and muscle growth (Tamate et al., 1962; Kristensen et al., 2007). With the development of high-throughput RNA sequencing (RNA-seq), it is allowed to investigate gene expression of certain tissue as integrity. In beef cattle, RNA-seq has been utilized to investigate the transcriptome of liver, longissimus dorsi muscle, and adipose tissues (Cesar et al., 2015; Huang et al., 2017; Kong et al., 2020). Of these, the developmental transcriptome of ruminal tissue in beef cattle was limitedly reported. Identifying genes expression during ruminal tissue development represents an important step toward understanding the biological processes of rumen growth.

Here, we performed RNA-seq to examine transcriptome profiles and identify candidate genes from rumen tissue of three developing stages. We believed that these results could help us understand the molecular mechanism of rumen development and further illuminate the relation between transcriptome of rumen development and beef complex traits.

Materials and Methods

Phenotypes and Rumen Tissue Collection

A total of 13 Simmental half-sib individuals from three periods (birth, youth, and adult) were used in this study, including five calves at birth, five youth individuals (6 months old), and three adult cattle (18 months old). These cattle were raised under the same feeding strategies and conditions in Shayang Hanjiang cattle Co., Ltd. (Shayang County, Jingmen City, Hubei Province). After slaughtering, a 2-cm2 piece of rumen tissue was isolated immediately, rinsed with sterilized PBS buffer (pH = 6.8) and placed in a 2-ml tube. All samples were immediately frozen with liquid nitrogen for total RNA extraction.

Nucleic Acid Extraction, Sequencing, and Genotyping

Total RNA was extracted from rumen tissue using TRIzol reagent (Invitrogen, Life Technologies) according to the protocol of instruction. Total RNA samples were assessed for RNA integrity using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 System (Agilent Technologices, CA, USA). RNA concentration was assayed by Qubit, and purity was assayed by Nanophotometer Spectrophotometer (Theermo Fisher Scientific, MA, USA). RNA-seq library construction was performed only for samples with RNA integrity greater than 7. The cDNA library was constructed using the Illumina TruSeqTM RNA Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. RNA-seq was performed on Illumina NovaSeq 6000 platform using a pairing end strategy [read length 150 base pairs (bp)].

Sequencing Data Analysis

The clean data were mapped to the reference genome (Bos taurus ARS-UCD1.2) by the HISAT2 (v2.2.1) (Kim et al., 2015). Effective reads aligned with gene regions were statistically calculated on the basis of genomic location information specified by bovine reference genome annotations (http://ftp.ensembl.org/pub/release-103/gtf/bos_taurus/). SAMtools (v1.9) (Li et al., 2009) was used to sort the BAM alignment files that were generated from HISAT2 by name. FeatureCounts (v2.0.1) (Liao et al., 2014) and Stringtie (v2.1.1) were used to estimate read counts and normalize reads as transcripts per kilobase of exon model per million mapped reads (TPM) for each sample (Pertea et al., 2015).

Differentially Expressed Genes Analysis

To investigate differentially expressed genes (DEGs) among three periods (birth vs. youth, birth vs. adult, and youth vs. adult), using DESeq2 (v1.30.1) to normalize the gene count data and calculate differential expression (Love et al., 2014). Genes with a Benjamini–Hochberg adjusted p-value < 0.05 were designated as differentially expressed.

Co-Expression Network Analysis

The DEGs from three comparisons were put together and removed redundant duplicate genes and using R package WGCNA to construct co-expression network (Langfelder and Horvath 2008). Co-expressed gene modules were detected by dynamic tree cutting method (Langfelder et al., 2008). Then, modules with highly correlated GS and MM values (p-value ≤ 0.05) and module–trait relationships with a correlation coefficient >0.5 were identified as stage-specific modules (Zhang et al., 2017).

Function Enrichment and PPI Analysis

To understand the function of genes in each stage-specific module, DAVID (https://david.ncifcrf.gov/) (Huang et al., 2007) and KOBAS (http://kobas.cbi.pku.edu.cn/kobas3/) (Xie et al., 2011) were used to GO enrichment analysis and KEGG pathway analysis, respectively. GO terms and pathways with a p-value less than 0.05 were defined as significantly enriched. Genes in each stage-specific module were calculated separately by STRING (https://string-db.org/) (Szklarczyk et al., 2015), obtained the protein–protein interaction (PPI) network, and imported into Cytoscape (v3.8.2). Cytoscape depends on the “Analyze Network” function to get the degree of a node, which means how many edges connect to it (Shannon et al., 2003).

The Enrichment Analysis of GWAS Signals

Our resource population included 1,478 Simmental beef cattle, which were born between 2008 and 2020 from in Wulagai, Inner Mongolia. All individuals were slaughtered at an average age of 20 months; four carcass traits [carcass weight (CW), lean meat weight (LMW), backfat thickness (BFT), and ribeye are (REA)], MS, and stomach weight (SW) traits were measured. Of these, SW is the total weight of the rumen, reticulum, and abomasum. All phenotypes were adjusted for the environmental fixed effects, including farm, year and sex, pre-fattening weight, and fattening days that were regarded as covariates. The DNA for each animal was isolated from blood samples and genotyped with Illumina BovineHD770kBeadchip (Illumina, San Diego, CA, USA). Before the statistical analysis, SNPs were filtered using PLINK (v1.90) (Purcell et al., 2007). SNPs were removed under strict criteria, and the standards are as follows: minor allele frequency (<0.01), missing genotypes (>0.05), and Hardy–Weinberg equilibrium (p < 10−6). Consequently, 1,432 individuals and 673,524 SNPs were remained (Xu et al., 2020). The mixed linear model was used for GWAS analysis by the GCTA (v1.93.0) (Yang et al., 2011; Yang et al., 2014) software, and the formula is as follows: where y is the pre-adjusted phenotypic value of the ith individual with the jth SNP; b is the allele substitution effect of SNP j; x is the jth SNP genotype of individual i and x is coded as 0, 1, and 2 for genotypes BB, Bb, and bb; g is the polygenetic effect of the ijth individual, g∼N (0, σa 2G), with σa 2 being the additive genetic variance and G is the additive genetic relationship matrix constructed using all SNPs. e is the residual effect, e ∼N (0, σe 2I), with σe 2 being the residual variance and I is the identity matrix.

In previous reported, rumen microbiota may affect the MS, adjusted 12th rib fat thickness, and other traits (Reddy et al., 2017; Kim et al., 2020), and its fermentation product VFAs can promote muscle growth (Kristensen et al., 2007). To investigate whether DEGs and gene co-expression modules were enriched with GWAS signals of carcass and other traits, we applied a sum-based method for GWAS signals enrichment analyses (sumGSE, https://github.com/WentaoCai/GWAS_enrichment) (Cai et al., 2021), which used signals of all markers within a pre-defined list of DEGs and then calculated the following summary statistics for the DEGs:

where Tsum was the summary statistics for a tested gene group. mg is the number of SNPs, which is within the DEGs or 5-kb up- and downstream of DEGs, and was acquired in GWAS statistics, meaning the value of marker effect. We randomly shifted the observed SNPs set to the new positions and calculated their Tsum summary statistics. The permutation was repeated 1,0000 times. When formula observed the proportion less than randomly sampled summary statistics, one-tailed tests was applied to calculate the empirical p-value.

Results

Transcriptome Profiling of Rumen Tissue

A total of 310.8 million raw reads from 13 rumen tissues were generated by RNA-seq. After quality control, we obtained an average of 23.1 million clean reads ranged from 19.1 to 29.0 million reads. The mapping rate was about 96.60% (ranging from 95.41% to 97.51%) after aligning clean reads to the cattle reference genome (ARS-UCD1.2) in Supplementary Table S1. These findings indicated good data quality that were suitable for subsequent analysis. We observed an average of 12,124 genes (ranging from 11,660 to 13,087) that were expressed (TPM > 0 in at least seven samples) across 13 samples. The gene expression of three periods is shown in Supplementary Figure S1, and the samples between youth and adult were highly correlated, while we found obvious differences between birth and youth/adult stage (Supplementary Figure S2). Despite differences in sample characteristics, samples from the same group were clustered together on the basis of their gene expression profiles (Figure 1A).

FIGURE 1

(A) PCA of the identified genes, the red, green and blue dots represent samples of adult, birth and youth periods. (B) Volcano plot of differential genes, volcano plot for DEGs in rumen tissue comparing birth period and youth period. (C) Volcano plot for DEGs in rumen tissue comparing birth period and adult period. (D) Volcano plot for DEGs in rumen tissue comparing youth period and adult period, red and blue dots represent up/down-regulated DEGs, respectively. The gray dots represent not DEGs.

Top Genes Expressed in Three Periods

The top 20 expressed genes in the rumen tissue at the birth, youth, and adult stages are shown in Table 1 and Supplementary Table S2. The highest expressed gene in rumen of birth stage was IRS4, GINM1, PAQR5, S100A14, MDP1, HOXC4, TAF7L, RS1, FDPS, NTS, JUP, PHB, UBE2G1, SMS, HBS1L, CGAS, TPM3, TUT4, GPR153, and GLOD4, whereas COX1, PJA1, S100A5, S100A12, ILF2, COX3, EFNB1, PCYT1B, USP11, CYTB, ND6, OPHN1, EDA2R, SPAG1, NEO1, ATP6, LCTL, COX2, RPLP1, JRKL and KRTAP11-1, COX1, DDX3X, S100A12, SLC39A1, SNORA75, COX2, COX3, ND3, ATP6, MID1IP1, ATP6AP2, PIP4K2C, RPLP1, ATP8, CXHXorf38, NPR1, NOX5, DES, and CYTB were highest expressed gene in youth and adult stage, respectively. We found that seven top expressed genes were same between youth and adult, including COX1, S100A12, COX3, CYTB, ATP6, COX2, and RPLP1, whereas the top 20 expressed genes of rumen in birth stage were totally different with the other stages, which indicated that the gene expressed pattern of rumen was distinct at birth stage.

TABLE 1

Top 20 expressed genes in the rumen tissues at three periods.

Period	Top 20 expressed genes
Birth	IRS4, GINM1, PAQR5, S100A14, MDP1, HOXC4, TAF7L, RS1, FDPS, NTS, JUP, PHB, UBE2G1, SMS, HBS1L, CGAS, TPM3, TUT4, GPR153, GLOD4
Youth	COX1, PJA1, S100A5, S100A12, ILF2, COX3, EFNB1, PCYT1B, USP11, CYTB, ND6, OPHN1, EDA2R, SPAG1, NEO1, ATP6, LCTL, COX2, RPLP1, JRKL
Adult	KRTAP11-1, COX1, DDX3X, S100A12, SLC39A1, SNORA75, COX2, COX3, ND3, ATP6, MID1IP1, ATP6AP2, PIP4K2C, RPLP1, ATP8, CXHXorf38, NPR1, NOX5, DES, CYTB

Differentially Expressed Genes Across Three Periods

Next, differences in gene expression of rumen tissues between different periods were investigated. In the birth vs. youth comparison, 4,905 DEGs were identified, including 2,486 upregulated genes and 2,419 downregulated genes (Figure 1B and Supplementary Table S3). A total of 3,877 DEGs were identified in birth vs. adult comparison, including 1,991 upregulated genes and 1,886 downregulated genes (Figure 1C and Supplementary Table S4). However, there were fewer DEGs identified between youth and adult, and only 521 DEGs were identified, including 314 upregulated genes and 207 downregulated genes (Figure 1D and Supplementary Table S5). The hierarchical clustering heatmap of all DEGs is shown in Supplementary Figure S3.

Co-Expression Network Construction and Module Detection

To better understand the function of DEGs, we used WGCNA to explore the relationship between DEGs. DEGs with a low expression level (TPM < 0.05) in more than one sample in the same group were removed, obtaining 6048 DEGs for co-expression analysis (Supplementary Table S6). A scale-free network was constructed using blockwise module function, and 6048 genes were divided into 14 modules, the number genes in each module, ranging from 116 in dark red module to 1,421 in turquoise module (Supplementary Figure S4 and Supplementary Table S7).

Period-specific Module Identification

GS and MM of all genes in the module were calculated to investigate the period-specific modules during rumen development (Supplementary Figure S5). GS was defined as the correlation between the gene and the developing period. MM was defined as the correlation between module eigengene and the gene expression profile. A strong correlation between GS and MM (p ≤ 0.05) shows that genes highly associated with a trait are often the most important elements of the modules associated with that trait. Furthermore, we use ME to represent the gene expression level in each module. The correlations between ME and the period of differentiation were analyzed (Figure 2). Ultimately, we identified 12 period-specific modules (average module–trait relationship > 0.5 and p ≤ 0.05), among which the modules royal blue, cyan, magenta, black, and turquoise were positively correlated with birth period; the green yellow, green, yellow, lightgreen, blue, and dark red modules were positively correlated with youth period; and the brown ones were positively correlated with adult period. In contrast, the modules yellow, lightgreen, blue, and dark red were negatively correlated with birth period, and the modules royal blue and turquoise were negatively correlated with youth and adult period, repectively.

FIGURE 2

Correlation between modules and differentiation period. The color, ranging from blue through white to red, indicates negative to positive correlation.

Visualization Potential Genes and Function Enrichment Analysis

Potential genes were obtained through PPI network and then visualized by Cytoscape (Supplementary Figure S6). The significantly enriched GO terms and the top 10 pathways of each module are shown in Figures 3 and 4 and Supplementary Figures S7 and S8. Detailed information of the GO terms and pathways is shown in Supplementary Tables S8, S9. For rumen development, the important pathways identified were butanoate metabolism, synthesis and degradation of ketone bodies, cell cycle, Wnt signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, and TGF-beta signaling pathway. Multiple significant GO terms are related to cell division, mitochondrial matrix, mitochondrial inner membrane, histone acetyltransferase complex, histone acetyltransferase activity, regulation of multicellular organism growth, positive regulation of phosphatidylinositol 3-kinase signaling, and positive regulation of MAPK cascade. Combining the expression level of DEGs, GO, and pathway results allows us to suggest CCNB1, CCNB2, IGF1, IGF2, HMGCL, BDH1, ACAT1, HMGCS2, and CREBBP as the promising candidate genes for rumen development. Potential genes related to rumen development and their enriched pathways and GO terms are shown in Tables 2 and 3.

FIGURE 3

Gene ontology (GO) enrichment analysis of related differentially expressed genes (DEGs). (A) GO enrichment analysis of DEGs in green module. (B) GO enrichment analysis of DEGs in cyan module. (C) GO enrichment analysis of DEGs in darkred module. (D) GO enrichment analysis of DEGs in black module. (E) GO enrichment analysis of DEGs in brown module. (F) GO enrichment analysis of DEGs in blue module.

FIGURE 4

KEGG pathway analysis of related DEGs. (A) The top 10 of pathway enrichment for green module. (B) The top 10 of pathway enrichment for cyan module. (C) The top 10 of pathway enrichment for darkred module. (D) The top 10 of pathway enrichment for black module. (E) The top 10 of pathway enrichment for brown module. (F) The top 10 of pathway enrichment for blue module.

TABLE 2

Description of potential DEGs that associated with rumen development from different module.

Gene name	Gene full name	Module
IGF1	Insulin-like growth factor I	Cyan
IGF2	Insulin-like growth factor II	Cyan
BDH1	3-Hydroxybutyrate dehydrogenase 1	Darkred
HMGCL	Hydroxymethyl-3-methylglutaryl-CoA lyase	Darkred
CCNB1	Cyclin B1	Lightgreen
CCNB2	Cyclin B2	Lightgreen
ACAT1	Acetyl-CoA acyltransferase 1	Lightgreen
HMGCS2	3-Hydroxy-3-methylglutaryl-CoA synthase2	Lightgreen
CREBBP	CREB Binding Protein	Turquoise

TABLE 3

Various GO terms and pathways related to rumen development shown from different modules.

GO term/pathway	ID	Gene	Module
Cell division	GO: 0051301	CCNB1, CCNB2	Lightgreen
Mitochondrial matrix	GO: 0005759	BDH1	Darkred
Mitochondrial inner membrane	GO: 0005743	HMGCL, HMGCS2, ACAT1	Darkred, lightgreen
Histone acetyltransferase complex	GO:0000123	CREBBP	Turquoise
Histone acetyltransferase activity	GO:0004402	CREBBP	Turquoise
Positive regulation of MAPK cascade	GO: 0043410	IGF2	Cyan
Regulation of multicellular organism growth	GO: 0040014	IGF1	Cyan
Positive regulation of phosphatidylinositol 3-kinase signaling	GO: 0014068	IGF1	Cyan
Cell cycle	bta04110	CCNB1, CCNB2, CREBBP	Lightgreen
Butanoate metabolism	bta00650	HMGCL, BDH1, ACAT1, HMGCS2	Darkred, lightgreen
Wnt signaling pathway	bta04310	CREBBP	Turquoise
MAPK signaling pathway	bta04010	IGF1, IGF2	Cyan
Synthesis and degradation of ketone bodies	bta00072	HMGCL, BDH1, ACAT1, HMGCS2	Darkred, lightgreen
PI3K-Akt signaling pathway	bta04151	IGF1, IGF2	Cyan
TGF-beta signaling pathway	bta04350	CREBBP	Turquoise

Different Group and Module Genes in the Enrichment of Association Signals

To investigate whether DEGs and gene co-expression modules were enriched with GWAS signals of carcass and beef quality traits, we applied enrichment analysis for all DEGs and modules across six traits (Tables 4, 5). As shown in Supplementary Figure S9A, DEGs were significantly enriched with GWAS signals of CW and LMW. In birth–youth comparison, fold change 0∼1 (slightly) upregulated DEGs were significantly enriched with GWAS signals of MS (Figure 5A), and fold change >3 (dramatically) downregulated DEGs were significantly enriched with GWAS signals of SW (Figure 5B). In birth–adult comparison, fold change 1∼2 up/downregulated DEGs were significantly enriched with GWAS signals of CW, LMW, REA, and BFT (Figures 5C,D). In youth–adult comparison, DEGs were not significantly enriched with GWAS signals. Furthermore, as shown in Supplementary Figure S9B, DEGs of turquoise module were significantly enriched with GWAS signals of CW, LMW, and REA; DEGs of lightgreen module were significantly enriched with GWAS signals of CW; and DEGs of other modules were not significantly enriched with GWAS signals (Supplementary Figures S9C,D). Of these, GOLGB1, ACVR2A, and TWIST2 belong to turquoise module and slightly upregulated DEGs, which indicated that GOLGB1 and ACVR2A might affect CW traits and TWIST2 might affect REA trait. Furthermore, IGF2 belongs to dramatically downregulated DEGs and might affect SW trait.

TABLE 4

The GWAS enrichment for DEGs in different groups and modules.

Trait	Group/module	p-value
A
CW	birth_adult	0.0039
LMW	birth_adult	0.0066
B
MS	up_fold_01	0.0031
SW	down_fold_314	0.0110
C
LMW	up_fold_12	0.0042
CW	up_fold_12	0.0020
REA	up_fold_12	0.0002
BFT	down_fold_12	0.0471
D
CW	Turquoise	0.0232
REA	Turquoise	0.0322
LMW	Turquoise	0.0296
CW	Lightgreen	0.0231

A. The three periods are compared in pairs. B. In birth_youth comparison, DEGs were divided into upregulated and downregulated genes, which were further divided into up/down_fold_01, up/down_fold_12, up/down_fold_23, and up/down_fold_3x according to fold change. C. In birth_adult comparison, the grouping principle is the same as that of B. D. DEGs were segmented by modules.

TABLE 5

Statistical description of six traits in Chinese Simmental beef cattle.

Traits	N	Mean (SD)	Min	Max
CW	1,471	284.40 ± 54.51	162.60	448.12
MS	1,436	5.19 ± 0.82	3.00	6.00
SW	1,059	8.45 ± 2.01	4.72	15.71
BFT	839	3.11 ± 2.68	0.05	12.50
REA	1,331	87.28 ± 14.39	51.00	133.00
LMW	1,456	240.26 ± 46.53	126.00	383.00

CW, carcass weight; MS, marbling score; SW, stomach weight; BFT, backfat thickness; REA, ribeye area; LMW, lean meat weight.

FIGURE 5

The GWAS enrichment for DEGs in different groups. (A) The enrichment for DEGs of four groups in birth vs youth. (B) The enrichment for DEGs of four groups in birth vs youth. (C) The enrichment for DEGs of four groups in birth vs adult. (D) The enrichment for DEGs of four groups in birth vs adult. The line means P is equal to 0.05.

Discussion

In this study, we obtained a comprehensive landscape of transcriptome profiles across 13 rumen tissue samples during three different stages of growth. Importantly, we identified candidate genes and networks related to rumen development. We obtained top 20 genes expressed in the rumen tissue during three periods. Of these, COX1, COX2, COX3, CYTB, ATP6, RPLP1, and S100A12 have occurred in both youth and adult period. COX1, COX2, and COX3 are three subunits of cytochrome c oxidase that may play an important role in the production of vast majority of ATP molecules in mammalian cells (Čunátová et al., 2020). CYTB is involved in electron transport in the mitochondrial respiratory chain (Seddigh and Darabi 2018). ATP6 has transmembrane transport activity of hydrogen ions (Li et al., 2016). RPLP1 plays an important role in the elongation step of protein synthesis (Du et al., 2007). S100A12 is involved in cell cycle progression and differentiation (Fritz et al., 2010). These genes and its functions indicated that two periods are associated with energy and cell proliferation. Furthermore, the top 20 expression genes of birth period had no common genes with other periods, but S100A14 in birth period and S100A5 and S100A12 in youth or adult periods belong to the same gene family that functions in cell cycle progression and differentiation (Liu et al., 2008).

Rumen development is regulated by a various of factors. Previous studies indicated the genes of FABP7, ILK, PDGFɑ, HMGCS2, FABP3, and AKR1C1 involved in rumen development (Kato et al., 2016; Malmuthuge et al., 2019). Moreover, CREBBP, TTF2, TGFB1, and PPARɑ are capable of contributing to the development of rumen epithelium (Baldwin et al., 2012; Connor et al., 2013; Connor et al., 2014). In our study, we identified that CCNB1, CCNB2, IGF1, IGF2, HMGCL, BDH1, ACAT1, HMGCS2, and CREBBP are the potential genes involved in rumen development, some of which were reported previously. HMGCS2 plays an important role in ketogenesis in the rumen epithelium of sheep during development (Lane et al., 2002). Furthermore, HMGCS2 was identified as a downstream target of PPARα. Increased production of VFA induced by intake of solid feed during weaning might promote ketogenesis in rumen epithelial cells. In this progress, PPARα promoted papillary development by activating HMGCS2 (Connor et al., 2013). ACAT1, BDH1, HMGCL, and HMGCS2 were found involved in the pathways of synthesis and degradation of ketone bodies and butanoate metabolism, which were found as features of rumen wall development (Sun et al., 2021). Microbial fermentation of dietary carbohydrates in the rumen produces large amount of short-chain fatty acids (SCFAs) (Poulsen et al., 2012), such as butyrate, which are known to affect rumen development and stimulate rumen epithelial cells proliferation in vivo (Sakata and Tamate 1978). SCFA transport process involves many regulatory factors including insulin-like growth factor (IGF), sodium hydrogen exchangers (NHE), monocarboxylate transporters (MCTs), and epidermal growth factor (EGF), which are also involved in the regulation of rumen epithelial cell proliferation (Baldwin 1999; Yang et al., 2012; Benesch et al., 2014; Nishihara et al., 2020). Butyrate is mainly absorbed and metabolized in the rumen epithelium (Sehested et al., 1999), producing ketone body under the action of ACAT1, BDH1, HMGCS2, and HMGCL (Kostiuk et al., 2010). Of these, HMGCS2 is a rate-limiting mitochondrial enzyme in the ketogenic pathway (Lane et al., 2002) and catalyzes synthesis of 3-hydroxy-3-methylglutaryl-CoA (HMG-CoA), the central metabolite of rumen epithelial cells (Xiang et al., 2016).

Cell proliferation is an important part of rumen development and is affected by various aspects (Naeem et al., 2012). In our study, CCNB1 and CCNB2 were enriched in the cell division, cell cycle pathway. Previous study indicated that the CCNB1 expression promoted ruminal cell cycle progression in goat (Gui and Shen 2016). CCNB2 was also found playing important roles in the acceleration of cell cycle and rumen development (Cui et al., 2018; Sun et al., 2021). CREBBP functions as a transcriptional coactivator of RNA polymerase II–mediated transcription and plays an important role in cell growth (Kalkhoven 2004; Baldwin et al., 2012). Furthermore, butyrate has been shown to affect the rumen by regulating the transcription factor CREBBP (Baldwin et al., 2012). CREBBP was enriched in the cell cycle, TGF-beta signaling pathway, and Wnt signaling pathway, part of which were highly associated with cell proliferation, apoptosis, and differentiation (Huang and Huang 2005; Kabiri et al., 2018; Guo et al., 2021). Kim et al. (2013) reported that the activity of IGF-I or IGF-II can regulate epithelial cell proliferation and differentiation in some tissues. Bach et al. (1995) suggested that IGF-II stimulates the proliferation and differentiation of rat myoblasts. Besides, these pathways are highly correlated with cell proliferation, apoptosis, and differentiation (Sun et al., 2015; Wang et al., 2018). In our study, IGF-I was enriched in the GO terms of positive regulation of phosphatidylinositol 3-kinase signaling and regulation of multicellular organism growth, and IGF-II was enriched in the GO terms of positive regulation of MAPK cascade. Hence, IGF-I and IGF-II may affect the proliferation and differentiation of rumen cells.

Rumen and its microbial abundance are closely related to the growth traits of ruminants (Andersen et al., 2021). Large amounts of VFAs and microbial proteins were produced in rumen with the assistance of rumen microbes (Krause et al., 2020; Newbold and Ramos-Morales 2020; Chai et al., 2021), which are important sources of energy, fatty acids, and amino acids for ruminants (Bergman 1990; Pathak 2008; Li et al., 2013; Liu et al., 2018). Fatty acids are associated with the deposition of marbling and backfat (Smith and Johnson 2015; Ueda et al., 2019; Zhang et al., 2019). Amino acids are the building blocks of proteins that are necessary for body protein synthesis such as muscle growth or milk protein secretion (Abdoun et al., 2006; Wu 2014). Previous research reported that SORT1, ITGA6, and TMEM39B were associated with intramuscular fat deposition (Cesar et al., 2018; Hérault et al., 2018; Ueda et al., 2021). Chromosome locations of PAQR3 co-located with QTL associated with fat deposition and one missense variant likely affect intramuscular fat via the LNPEP gene (Pena et al., 2013; Derks et al., 2021). We integrated differential genes obtained from the transcriptome with GWAS data using sum-based marker-set test method, which have been shown to be more potent or at least equivalent to most commonly used marker-set test methods for polygenic traits (Sørensen et al., 2017; Cai et al., 2021). In birth vs. youth comparison, the GWAS signals of MS trait were significantly enriched in fold change 0∼1 upregulated genes, which implied that these slightly upregulated genes were related to beef quality. Interestingly, we found that GWAS signals of SW trait were significantly enriched in fold change >3 downregulated genes, indicating that these dramatically downregulated genes may be correlated with weight of stomach. Furthermore, SORT1, ITGA6, TMEM39B, PAQR3, and LNPEP belong to fold change 0∼1 upregulated genes, and IGF2 belongs to fold change >3 downregulated genes. We can make a hypothesis that SORT1, ITGA6, TMEM39B, PAQR3, and LNPEP may be related to MS trait and IGF2 related to SW, which needs to be verified in future study. In birth vs. adult comparison, GWAS signals for CW, LMW, REA, and BFT were significantly enriched in fold change 0∼1 up/downregulated genes (Figures 5C,D). This implied that these slightly upregulated genes may be associated with CW, LMW, REA, and BFT traits. Many genes were reported to be associated with CW in previous studies. Fu et al. (2020) reported that a novel 65-bp indel in the GOLGB1 gene is related to chicken growth and body weight. PLA2R1 was associated with fat deposition and body weight, and the haplotype of ACVR2A gene had a significant effect on CW (Gheyas et al., 2015; Bhattacharya et al., 2016). Moreover, BMP7 and NFIA are associated with hot CW and LMW traits in beef cattle (Wang et al., 2020). TWIST2 and ACSL1 genes may be related to REA trait (Carvalho et al., 2019; Silva-Vignato et al., 2019). Positive correlations were found between FBP1 and PCCA gene expressions with BFT (Fassah et al., 2018). In addition, PECR and ACAT1 genes may be related to BFT trait (Piórkowska et al., 2017; Silva-Vignato et al., 2019). Moreover, these above 11 genes included in fold change 0∼1 up/downregulated genes group implied that these genes may be involved in regulating economic traits in beef cattle. This need for further investigation. Besides, the differential genes were divided into modules for enrichment analysis. We found that GWAS signals for CW, LMW, and REA were enriched in turquoise module, and GWAS signals for CW were enriched in lightgreen module (Supplementary Figures S8B,C). These implied that turquoise module genes might be related to CW, LMW, and REA traits and lightgreen module genes might be related to CW trait. We also found that GOLGB1, ACVR2A, and TWIST2 belong to turquoise module, which was consistent with our above hypothesis that GOLGB1 and ACVR2A might affect CW traits and TWIST2 affect REA trait.

Conclusion

This study explored the genes of CCNB1, CCNB2, IGF1, IGF2, HMGCL, BDH1, ACAT1, HMGCS2, and CREBBP on rumen proliferation and development. On the basis of our results, we can make a hypothesis that SORT1, ITGA6, TMEM39B, PAQR3, and LNPEP may regulate MS trait; GOLGB1, PLA2R1, and ACVR2A may regulate CW trait; IGF2 may regulate SW trait; BMP7 and NFIA may regulate LMW trait; TWIST2 and ACSL1 may regulate REA trait; and FBP1, PCCA, PECR, and ACAT1 may regulate BFT trait.