The Rumen Microbiome and Its Metabolome Together With the Host Transcriptome Act to Regulate the Cold Season Adaptation Mechanism of Grazing Tibetan Sheep

Background ： Tibetan sheep are important ruminants on the Qinghai-Tibet Plateau. They can maintain a normal life and reproduce in harsh environments under extreme cold and low oxygen. However, the molecular and metabolic mechanisms underlying the adaptability of Tibetan sheep during the cold season are still unclear. Hence, we conducted a comprehensive analysis of rumen epithelial morphology, epithelial transcriptomics, microbiology and metabolomics in a Tibetan sheep model to understand the interaction between the rumen host and microbiota and their metabolites and to explore the potential regulatory mechanism of Tibetan sheep adaptability to the cold season of the plateau. Results: Morphological analysis showed that the ruminal muscle layer thickness and nipple width of Tibetan sheep increased significantly during the cold season ( P <0.05), and the thickness of the stratum corneum, stratum granulosa and stratum spinous of the rumen epithelium increased significantly ( P <0.05). Transcriptomics analysis showed that the differential genes were primarily enriched in the PPAR signaling pathway (ko03320), legionellosis (ko05134), phagosome (ko04145), arginine and proline metabolism (ko00330), and metabolism of xenobiotics by cytochrome P450 (ko00980). Unique differential metabolites were identified in cold season, such as cynaroside A, sanguisorbin B and tryptophyl-valine, which were mainly enriched in arachidonic acid metabolism, arachidonic acid metabolism and linolenic acid metabolism pathways, and had certain correlation with microorganisms. Integrated transcriptome-metabolome-microbiome analysis showed that epithelial gene- GSTM3 expression was upregulated in the metabolism of xenobiotics by the cytochrome P450 (KO00980) pathway during the cold season, leading to the downregulation of some harmful metabolites; TLR5 gene expression was upregulated and CD14 gene expression was downregulated in the legionellosis (KO05134) pathway during the cold season. A large number of metabolites, such as glucosidic acid and vitamin A, were produced in the steroid hormone biosynthesis (KO00140) and retinol metabolism (KO00830) pathways. Conclusion: This study comprehensively described the interaction mechanism between the rumen host and microbes and their metabolites in grazing Tibetan sheep during the cold season. Our data show that under the stimulation of the cold plateau environment, the morphological structure of the rumen epithelium of Tibetan sheep undergoes adaptive changes. Rumen epithelial genes, microbiota and metabolites act together in some key pathways related to cold season adaptation. and the dominant positively the metabolites. Further study found that the were positively correlated but negatively correlated with creatine, acid,


Introduction
Tibetan sheep are primitive sheep bred under special environmental conditions on the Qinghai-Tibet Plateau (high cold, low oxygen and strong ultraviolet rays). Traditional Tibetan sheep breeding is primarily based on natural grazing. Grazing Tibetan sheep are widely affected by plateau ecological factors such as cold and nutritional stress during the cold season. The alpine grassland area of the Qinghai-Tibet Plateau can be divided into warm and cold seasons, corresponding to the grassy and the withed stage of grassland vegetation change, respectively. The grass starts turning green in early May and yellow in early October. The grassy period (warm season) lasts for less than 5 months. During the year, Tibetan sheep primarily obtain nutrients by eating natural forages, while during the long withered season (cold season), Tibetan sheep can only obtain nutrients by eating withered grass. The nutritional supply of Tibetan sheep during the warm and cold seasons is severely unbalanced, and they are under more severe nutritional stress during the cold season, leading them to form a growth pattern of "fat loss in cold season and rejuvenation in warm season". Their production efficiency is low, which seriously limits the utilization of germplasm characteristics and industrial development of Tibetan sheep.
Food digestion and absorption is a key process of animal adaptive evolution. In addition to the role of the animal's own genome, the intestinal symbiotic microbiome also plays an important role in helping the host digest food, synthesizing nutrients that cannot be synthesized by the host itself, thus expanding the animal's metabolic reservoir [1][2]. The host and microbial genomes must coordinate their work and perform their respective functions to maintain their body health under various environmental conditions [3]. Understanding the factors that affect the composition of microorganisms will help reveal the biological mechanisms by which microorganisms change the progress of individual characterization and clarify the nature of host-microbe interactions, which may lead to new coping strategies and methods. In addition, some human-based cross-sectional studies have identified external and internal effects, such as diet, environment, drugs, and host immune and metabolic status, but these effects account for only 10%-20% of the microbial diversity [4][5], and the diversity of most microorganisms among individuals is unexplainable. Therefore, more attention has been given to the integration of the host genome, transcriptome, epigenome, metabolome and microbiome interactions.
An increasing number of studies have shown the importance of uncovering the relevant laws and mechanisms underlying host physiological characterization by sequencing the gut microbiome and metagenomics comprehensively and correlating these data with the host genome, transcriptome and metabolic profile [6]. Based on multi-omics studies, the rumen microbiome and its metabolome and host metabolome have been shown to promote the individualized production performance of dairy cows jointly [7]; a comprehensive transcript and microbiome analysis showed that nutritional intervention improved the rumen function of stunted yaks and promoted compensatory growth [8]. The interactions among the calf rumen microbiome, rumen epithelial transcriptome and microbial metabolites suggest that a highly active early microbiome regulates rumen development at the cellular level, and miRNAs may mediate these host-microbial interactions [9]. At present, there have been many reports revealing the microbiota composition and gene function of the digestive tracts of ruminants living on the Qinghai-Tibet Plateau [10], but there is still a lack of information on rumen microbiota and the interaction of metabolites with the host in Tibetan sheep. We measured the rumen flora and rumen fermentation parameters of the grazing Tibetan sheep in the warm (July) and cold (December) seasons and found that the rumen SCFAs content and microbial flora abundance of Tibetan sheep in the cold season were significantly higher than those in the warm season; In the cold season, Bacteroidetes increased significantly, Firmicutes decreased significantly, while the abundance of Prevotella_1 and Succiniclasticum related to the production of propionic acid decreased, the Ruminococcus_2 related to fiber degradation increased significantly, and functions related to energy metabolism were significantly enriched. The interaction analysis between rumen microbe-SCFAs-host genes indicated that the differential expression of SGLT1, Claudin-4 and ZO-1 genes in rumen epithelium may be related to the lack of nutrients in cold season in Tibetan sheep, and SCFAs driven by rumen microbes are controlled and absorbed by epithelial-related genes [11], and the interaction between microbial flora density and SCFAs may be related to its adaptability to warm and cold seasons.
Therefore, we hypothesized that the rumen microbiome and metabolome actively participate in Tibetan sheep plateau adaptability during the warm and cold seasons through interactions with the host transcriptome. We used integrated bioinformatics-based 16S rRNA sequencing, a metabolome determination of rumen microbes and rumen transcriptome (mRNA sequencing) to explore the interaction between the host and microbes and their regulatory role in plateau adaptation during the cold and warm seasons. In addition, our detailed understanding of rumen development (function, morphology and colonization) during the warm and cold seasons may provide a means to use its functions to improve the productivity and health of ruminants in the future and to meet global food production needs.

Experimental design and sample collection
The grazing Tibetan sheep in Zuogaimanma Township, Hezuo City, Gannan Tibetan Autonomous Prefecture, Gansu Province, were taken as the research object, with an altitude of 3,300m. Twelve Tibetan sheep ewes (1 years old) with similar weight (35.12±1.43 kg) and good health conditions were randomly selected from the same herder's flock, all of which were in the local traditional natural grazing management state without any supplementary feeding. Samples were collected in the warm season (July) of 2019 and the cold season (December) of 2019. Before grazing in the morning, the rumen fluid was collected by sheep using gastric tube rumen sampler. The collected rumen fluid was quickly frozen in liquid nitrogen tank and brought back to the laboratory for storage at -80℃ for subsequent 16S rRNA analysis and metabolite identification. After the death of jugular vein bleeding, the rumen ventral sac tissue (1cm 2 ) was collected, the contents were rinsed gently with normal saline, and fixed in 4% paraformaldehyde for morphological analysis; At the same time, a small rumen ventral sac was cut, the contents were washed off quickly by PBS, and then the epithelial tissue was separated with blunt scissors, which was quickly stored in liquid nitrogen for subsequent total RNA extraction.

Morphological characteristics
The rumen ventral sac tissue was fixed in 4% paraformaldehyde for 24h, and then dehydration, transparency, waxing, embedding, slicing and dyeing. Hematoxylin and eosin were used for staining. Hematoxylin stained the nucleus blue-purple, and eosin stained the cytoplasm pink. The thickness of rumen muscle layer, papillary height, papillary width, cuticle thickness, granular layer thickness, spinous layer thickness and basal layer thickness were determined using CaseViewer section analysis system.

RNA extraction, transcriptome analysis and RT-qPCR verification
Trizol reagent method (DP762-T1C) was used to extract total RNA from the rumen epithelial tissue of Tibetan sheep, Nanodrop2000 was used for concentration detection. Agient2100, LabChip GX (platinum, Model Platinum Elmer LabChip GX) for integrity testing. The cDNA Library was constructed by the VAHTS Universal V6 RNA-SEQ Library Prep Kit for Illumina® Kit (NR604-02), and the Kit procedure was strictly followed. The products were further purified using a VAHTSTM DNA Clean Beads kit (N411-03). An Illumina NovaseQ6000 (San Diego) was used for the machine sequencing of the constructed library. A bioinformatics analysis was performed on the Biomark cloud platform BMKCloud (www.biocloud.net), and clean data were obtained after data filtering from raw data. Clean data were sequenced with the specified reference genome Ovis_aries (Oar_rambouillet_v1.0. Ovis_aries) by HISAT [12], and mapped data were obtained. StringTie [13] was used to compare reads on the pair for assembly. FPKM [14] (fragments per kilobase of transcript per million fragments mapped) was used to measure the transcription or gene expression level. The DESeq2 [15] data analysis method was used to analyze differentially expressed genes (the fold change (FC) represents the ratio of expression levels between two samples (groups), the false discovery rate (FDR) was obtained by adjusting the p value of significance), using fold change>= 2 and FDR<0.01 as the screening criteria for differentially expressed genes; and the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses of differential genes were performed by GOseq [16]. To verify the reproducibility and reproducibility of the gene expression data, the RNA-seq method was used to select 8 DEGs randomly from the individual RNA samples originally extracted by RNA-seq for reverse transcription-quantitative PCR (RT-qPCR). See Table  S1 (Additional file 2) for the primer information for these DEGs.

Microbial diversity analysis and LS-MS/MS metabolic profile determination
Microbial DNA extraction, high-throughput sequencing, and bioinformatics analysis have been described in previous reports [11]. Microbial species at the genus level were selected for the correlation analysis of this study using16S rRNA sequencing results.
Twelve rumen samples were analyzed with a liquid chromatography-mass spectrometry platform. After the sample was thawed at room temperature, 100 µL of sample was weighed each time, 500 µL of extract containing internal standard (1000:2) (volume ratio of methanol to acetonitrile=1:1, internal standard concentration 2 mg/L) was added, vortex mixed for 30 seconds, and then an ice water bath was followed by ultrasound for 10 minutes, followed by standing at -20℃ for 1 hour, and centrifuging at 4℃ for 15 minutes (12000 rpm). Then, 500 μL of supernatant was dispensed into an EP tube, and the extract was dried in a vacuum concentrator. A 150 μL extract (acetonitrile water volume ratio: 1:1) was added to the dried metabolite for resolution, followed by 30 seconds of vortexing, 10 min of ice water bath ultrasound, and 15 min of centrifugation at 4°C (12000 rpm). Lastly, 120 µL of the supernatant was removed into a 2 mL injection bottle, and 10 µL of each sample was mixed to form a QC sample for computer testing. The LC/MS system for metabolomics analysis was composed of a Waters Acquity I-Class PLUS Ultra High Performance Liquid Tandem Waters Xevo G2-XS QToF High Resolution Mass Spectrometer, and the column used here was purchased from a Waters Acquity UPLC HSS T3 column (1.8 µm 2.1*100 mm). The samples were eluted using positive (ESI+) and negative (ESI --) mobile phases consisting of water and 5% acetonitrile, 0.1% formic acid as solvent A and acetonitrile and 0.1% formic acid as solvent B at a flow rate of 0.35 mL/min and 400 µl/min. The subsequent mobile phase (A:B) elution gradient was0-0.25 min 98%-2%, 10.0-13.0 min 2%-98%, and 13.1-15.0 min 98%-2%, followed by ion source temperature: 150℃ and desolvation temperature: 500℃. The flow rates of the backblow and desolvent gas were 50 L/h and 800 L/h, respectively. The original data collected by MassLynx V4.2 were used for peak extraction, peak alignment and other data processing operations by Progenesis QI software. Based on the Progenesis QI software, online databases such as METLIN and self-built databases of BMG were used for metabolite identification according to sample types. The mass number deviation of fragment ion recognition is less than 100 PPM. BMKCloud (www.biocloud.net) was used to conduct a subsequent bioinformatics analysis on the identified metabolites. The differential metabolites were screened by combining the differential multiple, P value of the t-test and VIP value of the OPLS-DA model, and the screening standard was FC>1, P value<0.05 and VIP>1, and KEGG functional annotation and enrichment analysis were performed for differential metabolites.

Data analysis
The independent sample T test in the IBM SPSS Statistics 25 software was used for the statistical analysis of the morphological data, and the difference was significant and statistically significant with P<0.05; Spearman correlation test was used to analyze the correlation between genus-level rumen microbes (Top20) and differential metabolites and differential genes.

Morphological analysis of rumen epithelium
As shown in Table 1, the thickness of the muscular layer of the rumen wall was significantly greater during the cold season than it was in the warm season (P<0.01), the nipple width of the rumen epithelial layer was significantly greater in the cold season than in the warm season (P<0.05), and the nipple height in the warm season was significantly greater than that in the cold season (P<0.05). In addition, the rumen epithelial layer is divided into the stratum corneum (a), granular layer (b), spinous layer (c), and basal layer (d) ( Figure 1B, E) from the outside to the inside. Figure 1 shows that the rumen epithelium is different during cold and warm seasons. The thicknesses of the stratum corneum, granular layer, and spinous layer were significantly higher in the cold season than in the warm season (P<0.01), and the thickness of the basal layer in the warm season was significantly higher than that in the cold season (P<0.01). Figure 1 shows that there are differences in the musculature morphology of the rumen wall in Tibetan sheep during the cold season, and the thickness of the musculature during the cold season is significantly higher than that in the warm season (P<0.01), and the amount of connective tissue in the muscle layer in the cold season is greater than that in the warm season; The stratum corneum in the rumen epithelial layer during the cold season shows signs of partial shedding ( Figure 1A), and the shedding degree is greater than that in the warm season; in addition, the gap between the rumen papilla during the cold season is larger than that in the warm season, and the rumen papilla falls off to a greater degree.

Rumen epithelial transcriptome (RNA-seq) analysis 3.2.1 Analysis of differentially expressed genes
In this study, transcriptome sequencing analysis was performed on the rumen epithelial tissue of Tibetan sheep during cold and warm seasons, and 101.42 Gb clean data were obtained. The clean data from each sample reached 15.56 Gb, and the percentage of Q30 bases was 94.79% and above. The correlation diagram shows that the correlation values within the warm group were all greater than 0.907, and the correlation values within the cold group were all greater than 0.927. PCA showed that there was a significant difference between the warm group and the cold group, and the repeatability within each group was good. During the cold and warm seasons, 11,728 and 11,441 expressed genes were detected, of which 11,088 were co-expressed genes ( Figure 1C). Using fold change ≥ 2 and FDR <0.01 as the screening criteria, a total of 505 differentially expressed genes were found during the cold and warm seasons, of which 229 were upregulated during the cold season and 276 were downregulated ( Figure 1D). Compared with the warm season, the 5 genes with the highest expression in the cold season were TCHH, LOC114118453, GSTM3,

Functional enrichment analysis of differentially expressed genes
The differentially expressed genes were annotated, and 420 and 327 differentially expressed genes were annotated in the GO and KEGG databases, respectively. In the GO database, annotated genes are divided into three categories: biological processes (BP), cell components (CC), and molecular functions (MF). As shown in Figure 3, GO functional classification found that most genes were enriched in cells (GO:0005623), cell part (GO:0044464), organelle (GO:0043226), membrane (GO:0016020), and membrane part (GO:0044425). Second, it was enriched in cellular process (GO:0009987), single-organism process (GO:0044699), biological regulation (GO:0065007), metabolic process (GO:0008152), and response to stimulus (GO:0050896) genes. There were fewer differentially expressed genes annotated in MF, which were concentrated in binding (GO:0005488), catalytic activity (GO:0003824), molecular transducer activity (GO:0060089), and receptor activity (GO:0004872). Additionally, the statistics of GO enrichment results show that the number of differential genes enriched in BP (1370 GO terms) is the largest, followed by MF (457 GO terms) and CC (278 GO terms). As shown in Figure S1 (Additional file 2), the bubble diagram of BP enrichment shows that differentially expressed genes are primarily enriched in oxidation-reduction process (GO:0055114), RNA-dependent DNA biosynthetic process (GO:0006278), positive regulation of gene expression (GO:0010628), and regulation of the heart rate by cardiac conduction (GO:0086091).
KEGG enrichment analysis (q value <0.05) ( Figure 4) found that the differentially expressed genes were primarily enriched in the PPAR signaling pathway (ko03320), legionellosis (ko05134), phagosome (ko04145), arginine and proline metabolism (ko00330), and the starch and sucrose metabolism (ko00500) pathways. Among them, the PPAR signaling pathway is a pathway related to lipid metabolism and gluconeogenesis; the phagosome and legionellosis pathways are primarily related to the body's immune stress. Upregulated genes were significantly enriched in the metabolism of xenobiotics by cytochrome P450 (ko00980), chemical carcinogenesis (ko05204), steroid hormone biosynthesis (ko00140), retinol metabolism (ko00830) and other pathways; downregulated genes were primarily enriched in the PPAR signaling pathway (ko03320), estrogen signaling pathway (ko04915), AMPK signaling pathway (ko04152) and other pathways.  Up-regulation gene enrichment pathway analysis. c Down-regulation gene enrichment pathway analysis

RT-qPCR verification
Eight genes were randomly selected for verification, of which 4 genes were up-regulated and 4 genes were down-regulated. As shown in Figure 5, the qPCR expression patterns of the selected genes are consistent with the RNA-Seq analysis results, indicating the reliability and accuracy of the RNA-Seq method used in this study.

Analysis of rumen metabolome and microbiome 3.3.1 Differential metabolite analysis
The analysis of metabolites in the rumen during cold and warm seasons revealed 3,221 metabolites (positive model) and 3,104 metabolites (negative model). Principal component analysis (PCA) found that there were specific differences in the rumen metabolites of Tibetan sheep during cold and warm seasons under two different ion modes (positive and negative ions). Orthogonal partial least squares discriminant analysis (OPLS-DA) found that the R2X, R2Y and Q2 values were close to 1, where Q2>0.9, further verifying the reliability of the OPLS-DA Model (Additional file 2: Fig S2). FC>1, P value<0.05 and VIP>1 were used as the screening criteria to perform differential metabolite analysis (Fig. 6). In the positive model, a total of 1,649 differential metabolites were identified, with 824 upregulated and 825 downregulated differential metabolites, of which 25 were unique differential metabolites during the cold season (Additional file 3), such as hexanal, cynaroside A, and rishitin. In the negative model, a total of 1,599 differential metabolites were identified, with 1,048 upregulated and 551 downregulated differential metabolites, of which 131 were unique differential metabolites (sanguisorbin B, undecanedioic acid, geranyl acetone, tryptophyl-Valine, etc.). Further functional annotation of the differential metabolites showed that 145 (positive model) and 140 (negative model) differential metabolites were annotated to the KEGG function, and the annotations were classified and analyzed based on the HMDB database (Figure 7). We found 145 positive ionization differential metabolites, which were primarily classified into fatty acids (23), followed by carboxylic acids and derivatives (19), glycerophospholipids (15), etc. The 140 negatively ionized metabolites were primarily classified into glycerophospholipids (20), followed by carboxylic acids and derivatives (18) and fatty acyls (13).

Functional annotation and enrichment analysis of differential metabolites
In the positive mode, the differential metabolites between warm and cold seasons were primarily enriched in arachidonic acid metabolism, alpha-linolenic acid metabolism, linoleic acid metabolism, bile secretion and the metabolism of xenobiotics by cytochrome P450. In the negative mode, the differential metabolites were primarily enriched in purine metabolism, glycine, serine and threonine metabolism, folate biosynthesis, tryptophan metabolism and protein digestion and absorption (Additional file 2: Fig S3). The cold season-upregulated metabolic pathways enriched a total of 12 pathways in the positive mode ( Figure 8A), among which arachidonic acid metabolism, linoleic acid metabolism, and alpha-linolenic acid metabolism were the primary enriched metabolic pathways, involving hexanal, lipoxin A4, alpha-linolenic acid, gamma-linolenic acid (P=0.001904814), leukotriene A4 and other metabolites. A total of 23 metabolic pathways upregulated in the cold season were enriched in negative mode (Figurer 8B), including histidine metabolism; arginine and proline metabolism; phenylalanine, tyrosine and tryptophan biosynthesis; biosynthesis of amino acids (P=0.000643301); glycine, serine and threonine metabolism were significantly enriched metabolic pathways; and the differential metabolites primarily included histidine, arginine, proline, glycine, serine and threonine. During the cold season, downregulated metabolic pathways were enriched in a total of 9 in the positive mode. (Figigure 8C) The significantly enriched metabolic pathway was the bile secretion (P<0.001) pathway, and the primary differential metabolites were ibuprofen, lamivudine, leukotriene C4 and bilirubin. There were 27 metabolic pathways that were enriched in the negative mode ( Figure 8D). The metabolic pathways that were significantly enriched were glyoxylate and dicarboxylate metabolism (P<0.001), carbon metabolism (P=0.00001), and glycine, serine and threonine metabolism (P=0.00004). The different metabolites were primarily L-serine, glyceric acid, carbon dioxide, and 4-hydroxy-2-oxoglutaric acid. Fig. 8 Functional enrichment analysis of differential metabolites. a Functional enrichment diagram of up-regulated differential metabolites under positive model. b Functional enrichment of up-regulated differential metabolites in a negative model. c Functional enrichment of down-regulated differential metabolites in a positive model. d Functional enrichment of down-regulated differential metabolites in a negative model.

Interaction analysis of rumen metabolome and microbiome
The results of previous studies on microbial diversity [11] and metabolome correlation analysis were combined to reveal the relationship between metabolites and microbial taxon OTUs to analyze the microbial population structure, physiological metabolism and genetic variation. As shown in Figure 9, the correlation heat map reveals that there is a specific correlation between different metabolites and microorganisms. This study screened microorganisms and metabolites with a correlation C value of greater than 0.99 and a CCP<0.01 and found that Christensenellaceae_R-7_group, Ruminiclostridium_9, Zoogloea and the metabolites avocadene 2-acetate, lysyl-asparagine, gamma-glutamyl cysteinylserine, and biocytin have extremely significant and strong correlations (P <0.01). In addition, restricted correspondence analysis was performed for differential metabolites and microflora, as shown in the figure 9B. Meta2 (lumochrome), meta5 (D-pyroglutamic acid), meta_6 (L-pipecolic acid), meta_8 (caprylic acid (octanoic acid)), Meta_13 (inosine), meta_24 (theasapogenol E) and other differential metabolites are closely related to microorganisms.
As shown in Figure 9C, a correlation analysis between the top 20 microorganisms at the genus level and the differential metabolites was further performed. Bilirubin was extremely significantly positively correlated with Prevotellaceae_UCG-001 (P<0.001) and extremely significantly negatively correlated with Pseudobutyrivibrio (P<0.001). Creatine had a very significant negative correlation with Rikenellaceae_RC9_gut_group (P<0.001) and a positive correlation with the difference in Butyrivibrio_2 (P<0.05). Uncultured_bacterium_o_WCHB1-41 was positively correlated with bilirubin and negatively correlated with histamine. From the correlation heat map, it is clear that the dominant bacteria in the rumen of the cold season are positively correlated with the upregulated metabolites in the cold season, and the dominant bacteria in the rumen of the warm season are positively correlated with the downregulated metabolites. Further study found that the dominant bacteria Rikenellaceae_RC9_gut_group and carbon dioxide, bilirubin, 4-hydroxy-2-oxoglutaric acid, and taurolithocholic acid 3-sulfate were positively correlated but negatively correlated with creatine, guanidinoacetic acid, histamine, PC(20:5(5Z,8Z,11Z,14Z,17Z)/24:0) and chorismate.

Fig. 9
Correlation analysis of rumen microbes and metabolites. a Correlation analysis between differential metabolites and OUT. b Restriction correspondence analysis between differential metabolites and microflora. c Correlation heat map between differential metabolites and microflora. *P < 0.05, **P < 0.01, ***P < 0.001

Rumen transcriptome-microbiome-metabolome joint analysis
As shown in the figure, there was a correlation between the rumen microbial flora and the host transcriptome (Additional file 2: Fig S4); for example, the Rikenellaceae-RC9-gut-group, Prevotallaceae-UCG-003, Butyrivibrio-2, Ruminococcaceae-NK4A214-group and other flora had a significant correlation with host gene expression (P<0.05). In the transcript of the rumen epithelium, the differential genes that were upregulated in the cold season were significantly enriched in metabolism of xenobiotics by cytochrome P450 (KO00980), and in metabolome studies, differential metabolites were also found to be enriched in this pathway, as shown in Figure  10, due to the metabolic process of cytochrome P450, the expression of gen-GSTM3 is significantly upregulated in the cold season, and the metabolites S-(2,2-dichloro-1-hydroxy)-ethyl glutathione and 2-(S-glutathionyl)-acetyl chloride was significantly downregulated. In addition, the legionellosis (ko05134) pathway is involved in the immune stress response. Under the influence of the external environment and microbial metabolites, the expression of the TLR5 gene is upregulated in the cold season, and the expression of the CD14 gene is downregulated, which leads to the occurrence of proinflammatory response chemoattraction. In the microbial metabolome, the metabolites upregulated during the cold season were primarily enriched in the arachidonic acid metabolism and histidine metabolism pathways. The pathway diagram shows that linoleic acid metabolism produces arachidonate, which further produces LTA4, 15(S)-HETE, and LXA4, all of which are cold season upregulations. In the leukotriene metabolism pathway, it is metabolized into LTF4, which is downregulated during the cold season. In the histidine metabolism pathway, it participates in the pentose phosphate pathway, in which the cold-season downregulated metabolite AICAR nucleotide further participates in purine metabolism. L-Histidinol histamine alcohol was also upregulated and involved in the production of histamine metabolites as well as the production of Anserine metabolites. In addition, we found that cold season-upregulated genes are also enriched in the steroid hormone biosynthesis (ko00140) and retinol metabolism (ko00830) pathways, resulting in the production of a large amount of glycosides and vitamin A and other metabolites (Additional file 2: Fig S5/S6).

Fig. 10
Analysis of differential genes and metabolite regulatory pathways in rumen epithelium during cold and warm seasons. Boxes represent genes, circles represent metabolites, red represents up-regulation and green represents down-regulation.

Discussion
At present, the mechanism of interaction between rumen hosts and microorganisms and their metabolites in Tibetan sheep during cold and warm seasons is very limited. In this study, we integrated rumen epithelial morphology, epithelial transcriptomics, microbiology and metabolomics analyses to explore the interaction between the rumen host and microbiota and their metabolites, thereby revealing the regulatory mechanism of Tibetan sheep adaptability to the plateau cold season. The rumen epithelium is composed of leaf-shaped papillae, which not only serve as an absorption structure but also serve as an epithelial barrier to prevent the invasion of rumen microorganisms or toxins [17]. Recent studies have indicated that the characteristics (size) of fiber particles and the effective fiber content play a key role in rumen muscularization and volume development [18]; the higher the NDF (such as hay) is, the better the development of the rumen muscle layer [19]. In this study, the muscular thickness of the rumen wall in the cold season was significantly higher than that in the warm season, which was consistent with previous research results and might be related to the fact that Tibetan sheep grazed a large amount of withered grass during the cold season. On the one hand, the rumen muscle layer is adapted to the stimulation of the dry withered grass, which leads to an increase in the thickness of the rumen muscle layer; on the other hand, ingesting a large amount of withered grass that is not easily degraded and digested, the body needs to increase the level of rumen movement to degrade and digest the withered grass, thereby promoting the development of the rumen muscle layer. Physical stimulation from hay in the diet promotes the growth of rumen papilla, and the width of the rumen papilla is significantly increased by straw with a higher NDF [20]. In this study, under the stimulation of cold-season hay, the width of the rumen nipple in the cold season was significantly larger than that in the warm season, and the height of the papilla in the warm season was larger than that in the cold season. The rumen epithelium is composed of the stratum corneum, granular layer, spinous layer and basal layer and has important physiological functions such as VFA metabolism, absorption, transportation and protection [21][22]. The stratum corneum has a small metabolic capacity and is an important protective barrier [23], and this study found that the thickness of the stratum corneum during the cold season is significantly higher than that in the warm season. This finding may be due to the rough quality and high fiber content of the cold season. The rumen is strongly contracted and relaxed to complete mechanical digestion, which promotes an increase in the thickness of the stratum corneum, thereby forming a strong epithelial protective barrier. The basal layer is the most important rumen energy metabolism layer, and basal cells contribute to metabolic properties such as ketone production [24]. Granular layer cells establish tight gap junctions to maintain the integrity of the metabolite concentration gradient in the whole rumen wall [25]. The results of this study indicated that the thickness of the granular layer during the cold season was significantly greater than that in the warm season, while the basal layer appeared to be more compact in the cold season. These special morphological structures provided the possibility for Tibetan sheep to adapt to the cold season. To understand the molecular mechanism behind this morphological difference, rumen epithelial tissue was further analyzed by RNA-seq, and it was planned to explore its adaptive regulation mechanism during the cold season at the transcriptome level.
RNA-seq is an effective, efficient and widely used method for transcription status analysis and functional gene discovery [26][27][28]. In this study, the rumen transcriptome analysis of cold and warm seasons uncovered 505 differentially expressed genes, of which 229 were upregulated in the cold season and 276 were downregulated. TCHH was highly expressed during the cold season among the upregulated genes. Studies have found that TCHH and TCHHL2 have homology, and TCHHL2 plays a role in the cross-linking of keratin on the rumen surface [29]. In the abovementioned morphological studies, it was also found that the thickness of the stratum corneum increased during the cold season, which may be caused by the high expression of the TCHH gene. In addition, studies have found that GSTM3 has the effects of detoxification and the removal of reactive oxygen species (ROS) [30][31], and under the harsh environment of the cold season, to ensure the normal function of the rumen; it exhibits high GSTM3 gene expression. The study found that DECR1 encodes the rate-limiting enzyme for the oxidation of polyunsaturated fatty acids (PUFAs), which emphasizes the importance of DECR1 in the development of treatment resistance [32]. In this study, the DECR1 gene was highly expressed in the rumen during cold seasons. Based on the above analysis, these genes that were highly expressed in the cold season provide the possibility for Tibetan sheep to adapt to the cold season. In addition, we screened an important cold season-upregulated gene, LOC101113965. Studies have found that LOC101113965 (ubiquitin) plays an important role in regulating protein function. The dysregulation of the ubiquitin system leads to the occurrence of many diseases [33]. Ubiquitin proteins can also mark transmembrane proteins and participate in the transport of protein vesicles. The high expression in the rumen of the cold season regulates the protein in the rumen, thereby maintaining the normal function of the rumen. The further functional annotation analysis of its differential genes was performed, and in the GO database, most of the genes were enriched in processes related to cell growth and organelle growth. Lin et al. [34] also found that the differentially expressed genes of the rumen epithelium were primarily enriched in pathways related to cell growth, which is consistent with our research results. Under the influence of the lack of nutrients during the cold season, the growth and development of rumen epithelial cells are also affected to some extent, and the papilla morphology of the rumen epithelial layer has undergone adaptive changes. In addition, we also found that differentially expressed genes were enriched in some processes related to biological regulation and metabolism and in response to some external environmental stimuli. We speculate that under the stimulation of the cold season environment, some biological processes and metabolic functions of the rumen epithelium will change to make corresponding physiological and behavioral responses to adapt to the special cold season environment. In addition, the KEGG functional enrichment analysis of differential genes indicated that they were primarily enriched in pathways related to lipid metabolism, gluconeogenesis, and body immunity. Studies have suggested that the innate immune function of the rumen epithelium may be regulated differently during the transformation of the dietary structure [35], and changes in diet lead to changes in immune function, which is consistent with the results of this study. In the absence of nutrients in the cold season, the immune function of the rumen epithelium changes to a certain extent. Among them, the cold-season-upregulated genes are primarily enriched in the metabolism of xenobiotics by the cytochrome P450 (ko00980) pathway. In this pathway, we found that CYP1A1, GSTM and other genes are highly expressed, which leads to a higher ability of animals to eliminate rumen toxins [35]. These genes are reportedly involved in the metabolism of specific rumen toxins, such as ethanol and a series of exogenous substances [36][37], indicating that under the influence of the harsh cold season environment, the rumen organs of Tibetan sheep have a greater ability to remove toxins, thereby maintaining normal physiological functions. Second, upregulated genes are enriched in steroid hormone biosynthesis (ko00140) and retinol metabolism (ko00830). Studies have shown that steroid hormones mediate various important developmental and physiological functions in different organs, and a lack of steroids can cause adverse effects [38].
Adrenocortical hormone is a type of steroid hormone. It regulates glucose metabolism, inhibits the oxidation of sugar, raises blood sugar, and promotes the conversion of protein into sugar. Therefore, with the lack of nutrition during the cold season, the synthesis of steroid hormones ensures the body's normal sugar metabolism balance. Retinol (vitamin A) in retinol metabolism (ko00830) is a fat-soluble nutrient and an essential micronutrient. Vitamin A can increase incision inflammation and angiogenesis, repair collagen synthesis, and promote epithelial cell differentiation [39] to enhance intestinal resistance in the harsh environment of the cold season.
There is a strong symbiotic relationship between the animal intestinal microbiota and the host. The microbiota interacts with a variety of physiological functions in the host through its metabolites [40]. Therefore, we further used GC-MS to analyze the metabolic functions of the rumen microbiota in the cold and warm seasons and found that there were some differences in the rumen metabolites in the cold and warm seasons under two different ion modes (positive and negative ions). In the positive/negative mode, 25 and 131 different metabolites unique to the cold season were found, among which cynaroside A, sanguisorbin B, tryptophyl-valine and other metabolites were present in higher concentrations. Cynaroside is a flavonoid compound that is commonly found in honeysuckle [41]. Recent studies have shown that cynaroside is absorbed in the gastrointestinal tract, has an antioxidant effect, inhibits the oxidation of lipids and proteins [42][43][44][45], reduces cholesterol in atherosclerosis, and enhances capillary relaxation. In recent years, plant-derived bioactive compounds, such as saponins, have been considered potential alternatives to traditional antibiotics used as growth promoters in ruminant production [46]. Saponins have the ability to regulate rumen microbial communities and rumen metabolites. Research has found that elms have anti-inflammatory effects [47]. In addition, tryptophyl-valine is unique to the cold season. Studies have shown that valine is extensively metabolized in the rumen as a potential ketogenic and sugar-generating substance [48] and promotes growth [49]. We speculate that these unique metabolites play an important role in cold season adaptation. Further functional annotation classification of these differential metabolites revealed that both positive and negative ion compounds were classified into fatty acids and some of their derivatives, indicating that the differential metabolites in the cold and warm seasons were primarily fatty acid-related metabolites that participate in the body energy process. Among them, the positive ion metabolites that were upregulated during the cold season were primarily enriched in arachidonic acid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism and other pathways. Caldari-Torres et al. [50] found that the lipids of natural ruminant forages contain a large amount of essential linolenic acid, and linoleic acid can reduce cardiovascular disease [51]. Arachidonic acid is a biologically active substance of many cyclic eicosanic acid derivatives, and it is a ubiquitous component in mammalian cells. It is primarily found in glycerolipids or glycerophospholipids [52][53]. It is not only a tetra-unsaturated fatty acid that is important for normal cell membrane fluidity but also plays other important biochemical roles, including as direct precursors of bioactive lipid mediators, such as prostaglandins and leukotrienes [54]. These metabolic pathways were significantly upregulated in the cold season, suggesting that Tibetan sheep used essential fatty acids such as arachidonic acid, linoleic acid and linolenic acid to catabolize and adapt to nutritional stress and harsh environmental conditions in the cold season and to cope with nutritional stress and a low-temperature environment in the cold season. The anionic metabolites that were upregulated during the cold season were primarily enriched in some amino acid-related metabolic pathways. The amino acids in the rumen are the key precursors for protein and peptide synthesis and are primarily derived from dietary protein and trace protein [55], and rumen microbial proteins produced by rumen microorganisms can meet 90% of the amino acids that reach the small intestine [56]. Among them, arginine and proline are involved in RNA synthesis and protein glycosylation, which are necessary for cell function [57]. Proline is the main amino acid that maintains cell structure and function, and it is also an important regulator of cell metabolism and physiology. In addition, proline plays an important role in the synthesis and structure of proteins, metabolism and nutrition as well as the antioxidant response in trauma and immune responses [58]. In addition, the cationic metabolites that were downregulated during the cold season were primarily enriched in the bile secretion pathway. The differential metabolite bilirubin is the product of bile metabolism. Bile protects the body from intestinal infections by secreting immunoglobulin A (IgA) and inflammatory cytokines and stimulating the innate immune system of the intestines [59]. In the harsh environment of the cold season, the secretion of bile enhances the immune function of the rumen epithelium. According to research findings, the gut microbiota interacts with a variety of physiological functions in the host body through its metabolites [40]. Therefore, this study combined the results of previous microbial studies, provided a correlation analysis with the metabolome, and found that there is a correlation between different metabolites and microbes. The correlation analysis shows that Rikenellaceae_RC9_gut_group and Christensenellaceae_R-7 are associated with creatine and guanidinoacetic acid, respectively. It has a significant negative correlation with histamine and other metabolites. Studies have found that an increase in the abundance of the genus Christensenellaceae_R-7 enhances the catabolism of arginine [60]. Citrulline, a precursor for the synthesis of arginine, is positively related to Christensenellaceae_R-7 [60]. Arginine can synthesize guanidinoacetic acid and creatine [61]. This result further illustrates that the abundance of Christensenellaceae_R-7 is negatively correlated with the content of guanidinoacetic acid and creatine (p<0.05). The addition of saponins to the diet resulted in a decrease in the abundance of Prevotella_1 and Christensenellaceae_R-7, while alfalfa contained a large amount of saponins [60,62], further verifying that the abundance of Christensenellaceae_R-7 during the cold season increased significantly.
Based on the above analysis, we found that there is a correlation between the rumen epithelial transcriptome, microbiome and metabolome. To this end, this study used a combined analysis of the transcriptome-metabolome-microbiome to explore the molecular mechanism underlying the interaction between the host transcriptome and microbes and their metabolites during the rumen development of Tibetan sheep in cold seasons, thus revealing the special adaptability of Tibetan sheep to the cold season. The correlation analysis found that in the metabolism of xenobiotics by the cytochrome P450 (KO00980) pathway, the expression of gen-GSTM3 was significantly upregulated in the cold season, resulting in the metabolites S-(2,2-dichloro-1-hydroxy)-ethyl glutathione and 2-(S-glutathionyl)-acetyl chloride being significantly downregulated. The cytochrome P450 metabolic pathway plays a key role in the detoxification, cell metabolism and homeostasis of exogenous drugs [63]. Studies have found that GSTM3 belongs to the exogenous detoxification phase II enzyme family, and the detoxification of carcinogens is related to the metabolism of exogenous electrophilic substances. In addition, the legionellosis (ko05134) pathway is involved in the response to immune stress. Under the influence of the external environment and microbial metabolites, the expression of the TLR5 gene is upregulated in the cold season, and the expression of the CD14 gene is downregulated, which leads to the occurrence of proinflammatory response chemoattraction. TLR5 binds to bacterial flagellin to trigger the innate immune response to invading pathogens [64]. The human monocyte differentiation antigen CD14 is a pattern recognition receptor (PRR) that can enhance the innate immune response, and it has been confirmed that CD14 is a TLR co-receptor for detecting pathogen-related molecular patterns [65]. The enrichment of the legionellosis (ko05134) pathway in this study indicates that the innate immune mechanism of Tibetan sheep during the cold season has been improved. In the microbial metabolome, the metabolites upregulated in the cold season were primarily enriched in the arachidonic acid metabolism and histidine metabolism pathways. The pathway diagram shows that linoleic acid metabolism produces arachidonate arachidonic acid, which further produces LTA4, 15(S)-HETE, and LXA4, all of which are upregulated during the cold season. Lipoxin A4 (LXA4), one of the earliest endogenous lipid mediators, can inhibit the accumulation of neutrophils and inhibit the occurrence of inflammation [66]. In the histidine metabolism pathway, it participates in the pentose phosphate pathway, in which the cold-season downregulated metabolite AICAR nucleotide further participates in purine metabolism. L-Histidinol was also upregulated and involved in the production of histamine metabolites. Many tissues, especially the mast cells of the skin, lung and intestinal mucosa, contain large amounts of histamine. When tissues are damaged or inflammation and allergic reactions occur, histamine can be released, and vasodilation and vascular permeability can be generated through the H1 and H2 receptors [67]. These enrichment pathways could provide for the adaptability of Tibetan sheep to the cold season. In addition, we found that the cold season-upregulated genes were also enriched in the steroid hormone biosynthesis (ko00140) and retinol metabolism (ko00830) pathways, producing a large amount of glycosides and vitamin A and other metabolites (accessories), providing cold season Tibetan sheep with some energy requirements.

Conclusion
In summary, we integrated rumen epithelial morphology, epithelial transcriptomics, microbiology and metabolomics to analyze the interaction between the rumen host and microbiota and their metabolites, and the regulatory mechanism of Tibetan sheep adaptability in the cold season was revealed. Under the stimulus of the high-altitude cold season environment, the transcription status of the rumen epithelium and the rumen microbes and their metabolites of Tibetan sheep have undergone adaptive changes, resulting in changes in the morphological structure of the rumen epithelium. A transcriptomics analysis of the rumen epithelium indicated that differentially expressed genes were enriched in pathways related to cell and organelle growth in the rumen epithelium, which may lead to changes in rumen epithelial morphology. KEGG functional enrichment analysis indicated that differentially expressed genes were primarily enriched in pathways related to lipid metabolism, gluconeogenesis, and body immunity, thereby improving the immune function of Tibetan sheep during the cold season. The analysis of rumen microbes and their metabonomics showed that there is a correlation between differential microbes and metabolites, and cold-season differential microbes interact with a variety of physiological functions in the host through their metabolites. During the cold season, cynaroside A, sanguisorbin B, and tryptophyl-valine unique metabolites play a vital role in the antioxidation, anti-inflammatory, and growth processes of the rumen epithelium. Functional annotation analysis showed that these differential metabolites are all fatty acid-related metabolites, and they are primarily enriched in arachidonic acid metabolism, linoleic acid metabolism, alpha-linolenic acid metabolism, and amino acid metabolism-related pathways (arginine and proline metabolism), thus playing an important role in cell metabolism, physiological regulation and immune response. A comprehensive transcriptome-metabonomics-microbiome pathway analysis showed that during the metabolism of xenobiotics by the cytochrome P450 (KO00980) pathway, the expression of the epithelial gen GSTM3 was upregulated in the cold season, leading to the downregulation of some harmful metabolites. Therefore, the rumen organs of Tibetan sheep showed a higher ability to remove toxins physiologically. Under stimulation by the external environment and microbial metabolites, the expression of the TLR5 gene in the legionellosis (ko05134) pathway was upregulated in the cold season, and the expression of the CD14 gene was downregulated, which led to the improvement of the innate immune function of Tibetan sheep during the cold season. In addition to the steroid hormone biosynthesis (ko00140) and retinol metabolism (ko00830) pathways, a large amount of glycosides and vitamin A and other metabolites are produced, which provide for specific energy requirements for cold season Tibetan sheep. These key pathways provide insight into the molecular and metabolic mechanisms behind the cold season adaptability of Tibetan sheep.

Ethics Statement
All studies involving animal were carried out in accordance with the regulations for the Administration of Affairs Concerning Experimental Animal (Ministry of Science and Technology, China; revise in June 2004), and sample collection protocols were approved by the Livestock Care Committee of Gansu Agricultural University (Approval No. GAU-LC-2020-27).