Sixty 3-month-old post-weaned male Hu sheep with a healthy body condition, normal appetite, and similar body weight (24 kg–28 kg) were selected as the experimental animals. The animals were then randomly divided into three groups: the wheat straw group (WG), the alfalfa hay group (AG), and the peanut vine group (PG). Each group contained four replicates, and five sheep were included in one replicate. The nutritive information of the forage used in this study is provided in Table 1 . The feed composition and nutritive values of diet in each group are listed in Table 2 . The experiment was subdivided into a 7-day pre-trial stage and a 60-day trial stage. The diet nutritive levels were set in accordance with the Chinese industrial feeding standard of meat-producing sheep (NY/T 816-2004). All feedstuffs were mixed uniformly in a total mixed ration (TMR) machine and utilized instantly after mixing.

The daily management of the experimental animals was performed as usual, except that the conventional feed was gradually substituted by the diet we mentioned above in the pre-trial stage. During the 60-day trial stage, the experimental animals were fed twice a day (at 08:30 and 16:30) in an ad libitum manner, and they were also given free access to water. The animals’ daily feed intake (DFI) was recorded over the experimental period, and the ruminal fluid from one animal of each replicate was collected at the end of the feeding experiment, as was previously described (Zhu et al., 2022). In addition, the livestock status, including the animals’ appetite, water intake, and feces morphology, was monitored to ensure they were in healthy condition.

The average daily feed intake [ADFI, on a dry matter (DM) basis] was calculated by determining the difference between the weight of the feed provided to the animals and the weight of the residual feed for 1 day. The live weights of sheep were measured in the morning before providing feed and water at the first day and 60th day of the experiment. The total weight gain (kg) was calculated as the difference between the initial body weight (as measured on the first day) and the final weight (as measured on the 60th day). The average daily gain (ADG, g/d) was calculated by dividing the total weight gain (g) value by the duration of the experiment (60 days). The feed-to-gain (F/G) ratio was calculated by dividing the ADFI value by the ADG value.

On the 60th day of this experiment, four sheep from each group were selected, and 10 mL of blood was collected using a vacutainer from venous vessels in the neck. Subsequently, blood serum was obtained using a conventional method, that is, clotting at room temperature for 15 min, and centrifuging at 3,000 rpm for 15 min. After that, the serum was transferred into a 1.5-mL centrifuge tube and stored at a temperature of −20°C for later use. The serum levels of immunoglobin A (IgA), immunoglobin G (IgG), and immunoglobin M (IgM) were determined using a fully automatic biochemical analyzer (AU680, Beckman, CA, USA). The abundances of interleukin 1 beta (IL-1β; H002-1-2), interleukin 6 (IL-6; H007-1-1), interleukin 8 (IL-8; H008-1-1), tumor necrosis factor alpha (TNF-α; H052-1-2), interferon gamma (IFN-γ; H025-1-2), CD4 (H156-1-1), and CD8 (H157-1-1) were analyzed using the enzyme-linked immunosorbent assay (ELISA) method. All assay kits were purchased from Nanjing Jiancheng Bioengineering Institute (Nanjing, China). The experimental procedure was performed in accordance with the manufacturer’s guidelines.

On the 60th day of this study, one animal from each replicate was selected so that the slaughter performance of the sheep fed with the abovementioned three diet types could be assessed. The procedures of sample collection, slaughter performance examination, and meat quality determination (longissimus dorsi) were performed as previously reported (Jia et al., 2022).

After sacrifice, the compound stomach (including the rumen, reticulum, omasum, and abomasum) of the sheep was removed and washed to eliminate the attached chyme. Tissue samples measuring 2 cm × 2 cm in the middle of the ventral sac of the sheep rumen were harvested and stored in 10% neutral formalin solution for fixation. The measurement of the length, width, and surface area of ruminal papillae, and the thickness of submucosa, muscular, and mesothelium layers was carried out as reported before (Álvarez-Rodríguez et al., 2012).

The measurement of ruminal pH was performed instantly using a portable pH meter (Testo 206-pH2; Testo AG, Baden-Wurttemberg, Germany) after the sheep were sacrificed. Meanwhile, the ruminal fluid was filtered using four-layered gauze and stored at −20°C. The ruminal concentration of NH₃-N was measured using the classic phenol–sodium hypochlorite colorimetric method (Broderick and Kang, 1980). The concentration of SCFAs in the ruminal fluid was assessed using a high-performance liquid chromatography (HPLC)-based method (Isac et al., 1994).

The total DNA was extracted from rumen fluid samples using an E.Z.N.A.^® Stool DNA Kit (Omega Bio-Tek, Norcross, GA, USA). The V3 to V4 hypervariable region of the 16S rRNA gene was amplified using PCR with the following thermocycling protocol: 95°C for 2 min, followed by 30 cycles at 95°C for 20 s, 55°C for 30 s, and 72°C for 30 s, followed by a final extension stage at 72°C for 5 min. The primers used to amplify the sequences were as follows: V338F (5′-ACTCCTACGGGAGGCAGCAG-3′ and V806R (5′-GGACTACHVGGGTWTCTAAT-3′). The PCR reactions were performed in triplicate for 20-μL reactions consisting of 2 μL of 10 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 mM), 0.2 μL of FastPfu polymerase (TransGen Biotech, Beijing, China), 200 nM of each primer (Majorbio, Shanghai, China), and 10 ng of template DNA. The gel fragments of the correct size were excised and purified using an AxyPrep™ DNA gel extraction kit (Axygen, Union City, NJ, USA). The purified amplicons were equimolar pooled and paired-end sequenced (2 bp × 300 bp) on an Illumina MiSeq platform in accordance with standard protocols. All experimental procedures were performed in accordance with the manufacturer’s instructions. The raw reads were deposited into the National Center for Biotechnology Information (NCBI)’s Sequence Read Archive (SRA) database (accession number: PRJNA991257).

The raw 16S rRNA gene sequencing reads were demultiplexed, quality filtered using fastp (version 0.20.0) (Chen et al., 2018), and merged using FLASH (version 1.2.7) (Lozupone and Knight, 2005) in accordance with the following criteria: (i) the 300-bp reads were truncated at any site receiving an average quality score of < 20 over a 50-bp sliding window, the truncated reads shorter than 50 bp were discarded, and the reads containing ambiguous characters were also discarded; (ii) only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence; the maximum mismatch ratio of the overlap region was 0.2, and the reads that could not be assembled were discarded; and (iii) the samples were distinguished according to the barcode and primers, and the sequence direction was adjusted. The operational taxonomic units (OTUs) meeting the 97% similarity cutoff value were clustered using UPARSE version 7.1 (Edgar, 2013), and the chimeric sequences were identified and removed. The taxonomy of each representative sequence was analyzed using the RDP Classifier (version 2.11, http://sourceforge.net/projects/rdp-classifier/) against the 16S rRNA database (SILVA, version 138) using a confidence threshold of 0.7.

Pearson’s correlation analysis was performed using the online Majorbio platform (https://www.majorbio.com). The results were visualized using the R pheatmap package (version: 1.0.12). A p-value < 0.05 was set as significant for all analyses.

FastTree (version 2.1.3, http://www.microbesonline.org/fasttree/) was selected to construct the phylogenetic tree of the microbes at the genus level based on OTU data. The neighbor-joining method was adopted in the present study, and the result was visualized using R (version 3.3.1). The functional prediction of microbes was performed using the online tool PICRUSt2 (Douglas et al., 2020) on the Majorbio platform.

Reduced representation bisulfite sequencing (RRBS) was used to profile the methylome of the ruminal tissues in accordance with the standard protocol developed by Novogene Co., Ltd (Beijing, China). In brief, ruminal epithelial tissues were collected and rinsed twice to the remove the bacteria and adherent feedstuffs. Subsequently, the genomic DNA was extracted using the mammalian genomic DNA extraction kit (D0061, Beyotime, Beijing, China). A total of 1.5 μg of genomic DNA spiked with moderate lambda DNA was digested by MspI, after which it was subjected to end repair and adenylation. The cytosine-methylated barcodes were ligated to the DNA as per the manufacturer’s instructions. These DNA fragments were then treated with bisulfite using the EZ DNA Methylation-Gold™ Kit (Zymo Research). The sequencing libraries were constructed by Novogene Co., Ltd (Beijing, China). Subsequently, the paired-end sequencing of the samples was performed using the Illumina platform and 150-bp paired-end reads were generated (Illumina, CA, USA).

Subsequently, the obtained raw reads were assessed using FastQC (fastqc v0.11.5) and fastp (fastp 0.20.0) to acquire the filtered clean reads. The reads were then mapped to the sheep reference genome (ARS-UI Ramb v1.0) using Bismark software (version 0.16.3) (Krueger and Andrews, 2011). The reference genome was first transformed into a bisulfite-converted version (C-to-T and G-to-A converted) and then indexed using bowtie2 (Langmead and Salzberg, 2012). The clean reads were also transformed into fully bisulfite-converted versions (C-to-T and G-to-A converted) before being aligned to the similarly converted versions of the genome in a directional manner. The sequence reads that produced a unique best alignment from the two alignment processes (original top and bottom strand) were then compared with the normal genomic sequence, and the methylation state of all cytosine positions was inferred. The same reads that aligned to the same regions of genome were regarded as duplicates. The sequencing depth and coverage were summarized using the deduplicated reads. The results of the methylation extractor (Bismarck methylation extractor—no overlap) were transformed into the bigWig format for visualization. The sodium bisulfite non-conversion rate was calculated as the percentage of cytosines sequenced at the cytosine reference positions in the lambda genome. The methylated sites were identified by way of a binomial test using the methylated counts (mC), total counts (mC + umC), and the non-conversion rate (r). The sites with an FDR-corrected p-value < 0.05 were considered methylated sites. To calculate the methylation level of the sequence, we divided the sequence into multiple bins, with a bin size 10 kb. The sum of methylated and unmethylated read counts in each window was calculated. The methylation level (ML) for each window or C site shows the fraction of methylated Cs, and is defined as follows:

Differentially methylated regions (DMRs) were identified using the DSS software (Park and Wu, 2016). According to the distribution of DMRs throughout the genome, we defined the genes related to DMRs as the genes whose gene body region (from TSS to TES) or promoter region (2 kb upstream from the TSS) overlapped with the DMRs. The sequencing data were deposited in NCBI (accession number: PRJNA991642).

Gene Ontology (GO) enrichment analysis of genes related to DMRs was implemented using the GOseq R package (Young et al., 2010), in which the gene length bias was corrected. The GO terms with a p-value < 0.05 were considered significantly enriched by the differentially methylated promoter-related genes. We used KOBAS 3.0 software (Xie et al., 2011) to test the statistical enrichment of differentially methylated promoter-related genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

The data in the present study are presented as mean ± standard deviation (SD). One-way ANOVA was adopted to determine the statistical significance among the groups, and Duncan’s multiple comparison was used to discern the differences in the means among the three groups. A two-tailed Student’s t-test was used to compare the differences in the means between the two groups. A p-value < 0.05 was considered statistically significant.

To examine the impacts of roughage quality on the production performance of post-weaned Hu sheep, we assessed the growth performance, slaughter performance, and meat quality of animals fed low-quality (wheat straw) or high-quality (alfalfa hay and peanut vine) roughage-based TMR diets with nearly identical crude protein and digestive energy (DE) levels ( Table 1 ). The overall experimental design and associated downstream analysis are displayed in Figure 1 . As shown in Figure 2A , the ADFI of sheep in the wheat straw group (WG) was significantly lower than those of the animals in the alfalfa hay group (AG) and the peanut vine group (PG). Consequently, the ADG of animals in the WG was the lowest among the three groups, with sheep in the PG possessing the highest ADG ( Figure 2B ). Meanwhile, sheep in the AG and the PG displayed markedly lower levels of FCR than animals in the WG ( Figure 2C ). The above results suggest that animals fed high-quality roughage exhibit superior growth compared with those fed low-quality forage.

Subsequently, we demonstrated that the live weight of animals before slaughter was heavier in the PG than in the WG and AG ( Figure 2D ), with there being no significant difference existing between the WG and the AG (p > 0.05). Similarly, the dressing and net meat percentages of animals in the PG displayed elevated levels compared with those of animals in the WG ( Figures 2E, F ). At the same time, among the three groups, the bone-to-meat ratio of slaughtered animals was lowest in the PG ( Figure 2G ). Moreover, no significant differences were found in the internal organ index (e.g., heart), GR values, and lion eye muscle areas of slaughtered animals among the groups ( Supplementary Tables 1 and 2 ). These results indicated that the quality of the forage significantly affects the slaughter performance of post-weaned Hu sheep.

Finally, we tested the parameters pertaining to meat (longissimus dorsi) quality to dissect the impact of roughage grade on the eating characteristics of meat. As shown in Figure 2H , the drip loss of meat in sheep in the AG and PG was markedly lower than that of sheep in the WG (p < 0.05). The meat of sheep in the AG and PG was significantly more tender than that of sheep in the WG (measured by shearing force), although no statistical difference exists between the PG and WG ( Figure 2I ). In addition, among the three groups, the lightness (L*), redness (a*), and yellowness (b*) indices of the meat were highest in the PG, with slight differences observed between the WG and AG ( Figures 2J–L ). However, the pH_45min, cooked meat yield, and water-holding capacity of the meat were similar among the groups ( Supplementary Table 3 ). Collectively, these observations hinted that the roughage type exerts a substantial impact on the meat quality (i.e., eating quality and visual appearance) of Hu sheep.

To clarify the role of roughage quality on the health status of post-weaned Hu sheep, we determined the abundances of immunoglobulins, inflammation factors, and immune cells from peripheral blood. The levels of IgA, IgG, and IgM were significantly higher in the AG than in the PG and WG, with no statistical differences existing between the PG and WG ( Figures 3A–C ). The protein abundances of pro-inflammatory factors (i.e., TNF-α, IL-1β, IL-6, and IL-8) were significantly lower in the PG and AG than in the WG, with their levels in the AG being the minimum ( Figures 3D–G ). At the same time, the level of INF-γ in the AG is higher than those of WG and PG, even did not reach statistical significance standard (p > 0.05, Figure 3H ). Moreover, the numbers of CD4⁺ T lymphocytes ranked as follows: AG, PG, and WG ( Figure 3I ). However, the counts of CD8⁺ T lymphocytes were significantly higher in the AG than in the PG and WG ( Figure 3J ). The levels of CD4/CD8 were higher in the AG and PG than in the WG ( Figure 3K ). Notably, the quantities of CD4⁺ and CD8⁺ T lymphocytes, and CD4-to-CD8 ratio in the AG ranked first among the three groups. Taken together, these data clearly suggested that the fodder quality directly affects the partial key health status parameters of sheep, with the high-quality forage—alfalfa hay—exerting an optimally beneficial impact on the body condition of post-weaned sheep.

To further determine the influence of roughage quality on rumen fermentation and animal growth, we analyzed several parameters from ruminal fluids and rumen tissues. The pH of ruminal fluids in the three groups were nearly identical ( Supplementary Table 4 ). However, the concentration of NH₃-N significantly differed among the three groups, with that in WG sheep being the highest and that in AG sheep being the lowest ( Figure 4A ). The concentrations of SCFAs (acetate, propionate, butyrate, valerate, and total SCFAs) were higher in the PG than in the WG and AG ( Figures 4B–F ), whereas the ratios of A/P were lowest in the AG, and highest in the PG ( Figure 4G ). The ruminal papillae in the PG were dramatically longer than they were in the WG and AG ( Figures 4H, I ). However, the papillae width possessed an adverse trend compared with the papillae length within the groups ( Figure 4J ). The surface area of ruminal papillae was larger in the PG than in the WG and AG ( Figure 4K ). In addition, the rumen wall morphology was markedly influenced by forage quality, in that the submucosa of WG sheep was thicker than that of AG sheep and PG sheep ( Figure 4L ). The thicknesses of the muscular layer and mesothelium were greater in the AG than in the PG and WG ( Figures 4M, N ). Pearson’s correlation analysis results indicated that the pH of ruminal fluid was negatively associated with the papillae surface area ( Figure 4O ). The concentration of NH₃-N was negatively associated, to a high extent, with the mesothelium thickness, and, to a medium extent, with the thickness of the muscular layer. The concentrations of individual SCFAs and total SCFAs were generally positively linked to the parameters related to the rumen papillae and rumen wall. Notably, the concentration of butyrate showed an extremely high association with the thicknesses of the mesothelium and the muscular layer, and a slightly strong association with the thickness of the submucosa inside the rumen wall ( Figure 4O ). Generally, these findings implied that the roughage quality directly affects rumen fermentation, and potentially indirectly influences rumen development via SCFAs, especially butyrate.

To unearth the possible rumen microbiota-associated mechanisms underlying the differential production performances and rumen parameters among groups, we performed 16s rRNA sequencing of ruminal fluids, analyzed microbiota composition, and predicted the functional alterations associated with fluctuations in the microbiota. The rarefaction curves (indicated by the Sobs index at the OTU level) of all samples exhibited a steady pattern when the number of reads sequenced exceeded 10,000 ( Figure 5A ), thus indicating that the microbiome sequencing data were quantitatively adequate. Subsequently, we assessed the richness and diversity of microbiota in the three groups. The Aces index of the rumen fluid microbiome in the AG and PG was higher than of that in the WG, although it did not meet the statistical standard ( Figure 5B ). In contrast, the Simpson index values of the AG and PG were lower than that of the WG ( Figure 5C ). The above results hint that the roughage quality possibly alters the richness and diversity of ruminal fluid microbiota in post-weaned Hu sheep. We also performed β-diversity analysis (i.e., principal coordinate analysis, PCoA) to evaluate the general similarity of microbiota among groups. The obtained result hinted that the samples of WG are located in an individual area on the two-dimensional map, which clearly separates them from samples of the AG and PG ( Figure 5D ). The β-diversity analysis implies that the microbiota of the AG and PG share a higher similarity to that of the WG at the global level. The OTU data at the genus level are provided in Supplementary Table 5 .

Subsequently, we characterized the microbial composition of ruminal fluids collected from post-weaned Hu sheep at the genus level. The top ranked microbes (relative abundances > 1%) include: Prevotella 1, Rikenellaceae RC9 gut group, norank f F082, Ruminococcus 1, Ruminococcus 2, norank fMuribaculaceae, Ruminococcaceae NK4A214 group, Christensenellaceae R-7 group, Prevotellaceae UCG-001, and Succiniclasticum ( Figure 5E ). The majority of the microbes belonged to fibrolytic bacteria (e.g., Ruminococcus 1), amylolytic bacteria (e.g., Prevotella 1), secondary metabolite utilization-related bacteria (e.g., Succiniclasticum), and other microbes with undefined functions (e.g., Rikenellaceae RC9 gut group). We also identified the specific and shared microbes that existed in the three groups. As shown in Figure 5F , a large proportion of the microbes (189 genera) were common to all groups. A very small number of microbes were exclusively found in the WG (e.g., norank f norank o EMP-G18), AG (e.g., unclassified f Desulfovibrionaceae), and the PG (e.g., Amaricoccus). Furthermore, we defined the signature microbes of each group via intersecting the results from the ANOVA and Lefse analysis ( Supplementary Table 6 ). As shown in Figure 5G , seven microbes (including Succiniclasticum, Saccharofermentans, Family XIII AD3011 group, and Papillibacter) were filtered as the marker microbes for the WG. The Eubacterium ventriosum group, Erysipelotrichaceae UCG-004, and norank f Prevotellaceae were found as signature bacteria of the AG. Erysipelotrichaceae UCG-008 was screened as the marker species of the PG. The relative abundances of the first-ranked microorganism for each group are displayed in Figure 5H . We also performed Pearson’s correlation analysis of the relative abundances of microbes with the concentrations of SCFAs and other parameters ( Figure 5I ). The results suggested that the abundances of Ruminococcaceae UCG-014 are negatively correlated with the A/P ratio. The growth parameters (i.e., ADFI and ADG) showed negative relationships with the abundances of several microbes, including unclassified f Veillonellaceae, Veillonellaceae 001, Prevotellaceae UCG-003, Selenomonas 1, Succiniclasticum, Prevotellaceae UCG-001, norank f F082, and Rikenellaceae RC9 gut group. Several microbes (e.g., Rikenellaceae RC9 gut group) were among the top-ranked organisms of the ruminal fluid microbiota across all samples ( Figure 5E ). In addition, the level of Ruminococcaceae UCG-001 was positively associated with the levels of all SCFAs other than propionate, suggesting a potential pivotal role of the microbe in microbial decomposition of diet carbohydrates and production of partial SCFAs. Moreover, the level of NH₃-N exhibited a positive correlation with the norank f Bacteroidales RF16 group, suggesting that the microorganism participates in the microbial degradation of proteins, polypeptides, and amino acids. We also noticed that the abundance of Christensenellaceae R-7 group was negatively related to the content of propionate, and that the same relationship exists among Rikenellaceae RC9 gut group and the levels of butyrate and valerate. These analyses hinted that ruminal microbes are involved in the regulation of the rumen fermentation and animal growth, and that these relationships are influenced by forage quality.

Finally, we adopted PICRUSt2 (Douglas et al., 2020) to forecast the functional information of microbiota inside ruminal fluids. As shown in Figures 5J, K , several enzymes (3.1.3.18, phosphoglycolate phosphatase; 3.6.4.12, DNA helicase; and 6.3.1.2, glutamate–ammonia ligase) involved in nutrient degradation and synthesis, and five modules related to carbohydrate metabolism [e.g., M00001, glycolysis (Embden–Meyerhof pathway), glucose → pyruvate], were functionally enriched. These results implied that the forage quality changes the microbial functions related to carbohydrate metabolism. All the functional prediction results are provided in Supplementary Table 7 .

To understand the interrelationship of microbes in the rumen ecosystem and speculate their significant functions in rumen metabolism, we performed clustering, phylogenic, and correlation analyses using the OTU data of the microbes. As shown by the heat map in Figure 6A , the microbes with similar abundances across three groups clustered together, suggesting that these microbes possess similar functions or work synergistically to decompose polymers in the rumen. For example, Ruminococcus 1 (containing the main cellulolytic species Ruminococcus flavefaciens and Ruminococcus albus) and Ruminococcus 2 were gathered together, indicating that Ruminococcus 2 plays an important role in fiber degradation. Meanwhile, the succinate-producing bacteria (Prevotellaceae UCG-001) and the succinate-utilizing bacteria (Succiniclasticum) emerged in the same minor clade, substantiating the possible existence of a cross-feeding relationship between them. Moreover, the correlation analysis results of the bacterial abundances verified the above discoveries ( Figure 6B ). The abundances of microbes with similar functional characteristics or close metabolic interdependent relationships in the rumen were positively correlated. At the same time, the microorganisms with competitive relationships were negatively associated. For instance, the abundances of Prevotella 1 and Prevotella 7 showed a positive correlation relationship, whereas the abundance of a cellulolytic bacterium (Ruminococcus 1) displayed a negative association with an amylolytic microbe (Prevotellaceae UCG-001). The above results indicated that the abundance fluctuations of ruminal bacteria were consistent with their respective roles or complex metabolic relationships in maintaining the ecosystem of rumen.

The phylogenic analysis revealed that the top 20 ranked microbes are mainly divided into four clades ( Figure 6C ). Clade 1 is mainly composed of the bacteria belonging to the Ruminococcaceae family. Clade 2 comprises the microbes belonging to the Prevotellaceae family. The microorganisms in clade 3 belong to diverse families but are in the same order (Bacteroidales). A similar observation was made for the microorganisms in clade 4, except that all the microbes are in the Selenomonadales order. The above phylogenetic study hinted that the functions of microbes in the same clade are potentially relevant.

We also visualized the relationship of microbes at a global level through cooccurrence analysis and found that a complex interaction network exists. As shown in Figure 6D , the abundances of the microbes in group 1 (Ruminococcaceae family) were generally negatively associated with those of the microbes in group 3 (Prevotellaceae family) and those in Succiniclasticum. In addition, the levels of four microbes in group 2 were positively linked. However, the abundances of Rikenellaceae RC9 gut group and norank fF082 were negatively related to Prevotella 7. We also visualized the interactive map of microbes with rumen fermentation parameters and growth performance, and found that the microbes in group 2 are negatively correlated with AGD and ADFI ( Figure 6E ). The microorganisms in the groups contained amylolytic bacteria (e.g., Prevotellaceae UCG-001), intermediate product-utilizing bacteria (e.g., Succiniclasticum), and bacteria with unknown functions in the rumen (e.g., Rikenellaceae RC9 gut group). Moreover, the relationships of several bacterium with SCFAs and NH₃-N were displayed. The above results demonstrated that certain microbes are potential regulators of rumen fermentation and animal growth.

Based our previous observation that the rumen development is affected by fodder quality, we performed RRBS to explore whether or not the universal epigenetic status of rumen tissues is affected by that. We obtained highly qualified data from the rumen tissues of animals in the PG and WG ( Supplementary Table 8 ). Subsequently, we found that the methylation levels of the samples were highly correlated at a global level ( Figure 7A ). However, the proportions of methylated nucleotides exhibited an elevated trend in the gene body, upstream 2K, and downstream 2K regions in the PG compared with the WG ( Figure 7B ). Moreover, the methylation levels of differentially methylated regions (DMRs) in the PG were higher than in the WG ( Figure 7C ). Because the methylation status of gene promoter regions is closely related to gene expression via transcriptional regulation, we focused on the promoter segments of DMRs. A total of 164 genes with hypermethylated promoters and 117 genes with hypomethylated promoters were identified in the PG compared with the WG ( Figure 7D ; Supplementary Table 9 ). The KEGG pathway enrichment results hinted that the FoxO signaling pathway, aldosterone-regulated sodium reabsorption, T-cell receptor signaling pathway, carbohydrate digestion and absorption, VEGF signaling pathway, mTOR signaling pathway, and the sulfur relay system were overrepresented in the genes with hypermethylated promoters ( Figure 7E ). Moreover, the microRNAs in the cancer and mRNA surveillance pathways were found to be enriched genes with hypomethylated promoters. GO analysis suggested that the majority of genes with hypermethylated promoters are involved in organism development (e.g., anatomical structure morphogenesis), transcription factor functionality (e.g., DNA-binding transcription factor activity), and transepithelial water transport ( Figure 7F ; Supplementary Table). The enriched GO terms of the genes with hypomethylated promoters included cell junction, synapse, multicellular organism development, and system development. The detailed results of gene functional enrichment analysis are provided in Supplementary Table S10 . Furthermore, we demonstrated that the promoter regions of ACADL and ENSORAG0020014533 are hypermethylated and hypomethylated, respectively ( Figures 7G–J ). These results hinted that forage quality possibly determines the rumen development via altering the rumen methylation status and gene expression at an epigenetic level.