Histological analysis of the rumen
Results of H.E. staining of rumen tissue sections were presented in Figure 1A. The papillae length in the L group was significantly shorter than that in the H group (P < 0.05). Both papillae width and muscular layer of H group was significantly wider than in the L group (P < 0.05). No obvious difference was observed in the stratum corneum between the groups (Figure 1B, P > 0.05).
Antioxidant and immune response index assays
As shown in Figure 2A, the activities of T-AOC and SOD in rumen of H group were significantly higher than those of L group (P < 0.05). There was no significant difference in CAT, GSH-Px activity and MDA content in rumen fluid between the two groups (P > 0.05). Although the difference was not significant, there was an increased tendency of IgA, IgM, IgG, IL-6, and IL-1β in the H group than the L group (Figure 2B, P > 0.05).
RNA-Seq and Mapping
Table 1 showed the basic statistics for RNA-seq reads generated from the rumen of Tibetan sheep. A total of 390.63 million raw reads were obtained through sequencing in ten samples. After quality control, an average of 65.01 million clean reads was generated (64.04 million for H group and 65.97 million for L group). For all samples, the percent of clean reads was 99.85%. The Q30 values were between 94.01% and 96.75%. Average of 96.71% of the clean reads was mapped to the Ovis aries reference genome (H group with 97.52%, L group with 95.91%).
Analysis of Gene Expression
Venn diagram showed that there were 1706 and 893 differentially expressed genes (DEGs) belonging to this item in treatment and control DEG datasets, respectively, and they shared 14,458 DEGs as the intersection (Figure 3A). According to the gene expression values (FPKM), the principal component analysis (PCA) exhibited that H group and L group was separated into two different groups (Figure 3B). Hierarchical clustering analysis of the DEGs showed that the same group of DEGs was clustered together, indicating the accuracy and reliability of samples (Figure 3C). A total of 612 significant DEGs (based on |(fold change)| > 1.2, P-value < 0.05, and false discovery rate (FDR) < 0.05) were found in response to dietary CP level. Of these genes, 20 were up-regulated and 426 were down-regulated in the H group (Figure 3D).
Enrichment Analysis of DEGs
Gene ontology (GO) enrichment analyses were carried out using the database for annotation, visualization and integrated discovery software (https://david.ncifcrf.gov) (Figure 4A). A total of 573 significantly enriched GO terms were annotated within three major function groups: biological process (BP, 354), cellular component (CC, 95), and molecular function (MF, 124). For the down-regulated terms, the most significant GO categories observed were perikaryon, long-chain fatty-acyl-CoA catabolic process, and endoplasmic reticulum polarization. For the up-regulated terms, the most significant GO categories observed were potassium ion homeostasis, muscle contraction, and cell periphery.
Based on Kyoto Encyciopedia of Genes and Genome (KEGG) pathway enrichment analysis, 33 KEGG pathways were significantly (P < 0.05) enriched in the rumen tissue (Figure 4B). The top five pathways with the most representation of DEGs were retinol metabolism (13 DEGs), metabolic pathways (56 DEGs), ECM-receptor interaction (8 DEGs), focal adhesion (12 DEGs), and protein digestion and absorption (8 DEGs). Compared to the L group, 13 differentially nutrition metabolism-related signaling pathway (i. e., protein digestion and absorption, nicotinate and nicotinamide metabolism, glutathione metabolism, carbon metabolism, and fatty acid biosynthesis) were identified in H group.
Nine significant DEGs, identified from the RNA-seq data, were randomly selected for RT-qPCR validation, including two fatty acid biosynthesis-related genes (ACACB and ACSF3), seven nutrition metabolism-related genes (ATP1B1, IDH2, NT5E and NNT) and three muscle development-related genes (MYH11, KCNMA1 and MYL9). The RT-qPCR confirmed that the DEGs exhibited the same pattern of expression as observed with the RNA-seq (Figure 5), indicating our transcriptomic analysis was highly credible.
Correlation analysis
Correlation analysis (Figure 6) showed that T-AOC was positively correlated (P < 0.01) with ACACB (0.93), ACSF3 (0.94), ATP1B1 (0.94), and MYH11 (0.92), and IgG was positively correlated (P < 0.01) with ACSF3 (0.93). T-AOC was positively correlated (P < 0.01) with IDH2 (0.91), NT5E (0.86), KCNMA1 (0.91), MYL9 (0.87) were positively correlated (P < 0.05). IgG was positively correlated with ACACB (0.90), ATP1B1 (0.88), IDH2 (0.87), ATP1A2 (0.90), MYH11 (0.92), KCNMA1 (0.86) (P < 0.01). IL-6 was positively correlated (P < 0.05) with ACACB (0.83), ATP1B1 (0.87), IDH2 (0.84), ATP1A2 (0.83), MYH11 (0.95), KCNMA1 (0.85), MYL9 (0.84).