Enhanced propionate and butyrate metabolism in cecal microbiota contributes to cold-stress adaptation in sheep

Abstract

Background

During cold stress, gut microbes play crucial roles in orchestrating energy metabolism to enhance environmental adaptation. In sheep, hindgut microbes ferment carbohydrates to generate short-chain fatty acids (SCFAs) as an energy source. However, the mechanisms by which hindgut microbes and their metabolites interact with the host to facilitate adaptation to cold environments remain ambiguous. Herein, we simulated a winter environment (− 20 °C) and provided a rationed diet to compare the cold adaptation mechanisms between Hulunbuir and Hu sheep.

Results

Our findings show that cold exposure enhances SCFA metabolism in the sheep cecum. In Hu sheep, acetate, butyrate, and total SCFA concentrations increased, whereas in Hulunbuir sheep, propionate and butyrate concentrations increased, with a notable increase in total SCFAs. Notably, butyrate concentration was higher in Hulunbuir sheep than in Hu sheep under cold stress. Following cold exposure, the proinflammatory cytokine IL-1β levels increased in both breeds. In addition, Hu sheep showed increased IL-10, whereas Hulunbuir sheep exhibited elevated secretory IgA levels. The cecal microbiota responded differently, Hu sheep showed no notable changes in alpha and beta diversity, whereas Hulunbuir sheep exhibited considerable alterations. In Hu sheep, the abundance of fungi, specifically Blastocystis sp. subtype 4, decreased, and that of several Lachnospiraceae species (Roseburia hominis, Faecalicatena contorta, and Ruminococcus gnavus) involved in SCFA metabolism increased. Pathways related to carbohydrate metabolism, such as starch and sucrose metabolism, galactose metabolism, and pentose and glucuronate interconversions, were upregulated. In Hulunbuir sheep, the abundance of Treponema bryantii, Roseburia sp. 499, and Prevotella copri increased, with upregulation in pathways related to amino acid metabolism and energy metabolism. Cold exposure increased node connectivity within the symbiotic networks of both breeds, with increased network vulnerability in Hu sheep. Following cold exposure, the microbial community of Hulunbuir sheep showed a decrease in the influence of stochastic processes on community assembly, with a corresponding increase in the role of environmental selection. Conversely, no such shift was evident in Hu sheep. Further transcriptomic analysis revealed distinct regulatory mechanisms between breeds. In Hu sheep, protein synthesis, energy metabolism, and thermogenesis pathways were substantially upregulated. By contrast, Hulunbuir sheep showed considerable upregulation of immune pathways and energy conservation through reduced ribosome synthesis. Correlation analysis indicated that butyrate holds a central position in both networks, with Hulunbuir sheep exhibiting a more complex and tightly regulated network involving SCFAs, microbiota, microbial functions, and transcriptomes. Partial least squares path modeling revealed that cold exposure substantially altered the cecal microbiota and transcriptomes of Hulunbuir sheep, affecting SCFAs and cytokines.

Conclusions

The findings of this study suggest that under cold exposure, Hu sheep enhance acetate fermentation and rely on tissue thermogenesis for adaptation. By contrast, Hulunbuir sheep exhibit changes in microbial diversity and function, leading to increased propionate and butyrate metabolism. This may promote physiological energy conservation and innate immune defense, balancing heat loss and enhancing cold adaptation.

Introduction

Climate change poses a notable challenge to sustainable livestock development. Livestock farming must endure harsh winter conditions, including rain, snow, strong winds, and low temperatures. Extreme cold weather events threaten livestock farming because cold stress weakens animals’ resilience to natural disasters, affecting growth and health. This results in slow growth, weight loss, low lambing rates, and high mortality [1,2,3,4,5,6]. In addition, cold stress can lead to various diseases, such as respiratory infections [7, 8] and intestinal inflammations [9,10,11,12], severely affecting the economic benefits of farming.

Under cold stress, animals mobilize tissues and organs through the neuroendocrine system to increase heat production and reduce heat loss, thereby maintaining homeostasis. Prolonged cold exposure shifts heat production from shivering thermogenesis to nonshivering thermogenesis. The gut microbiota plays a crucial role in regulating the host’s metabolic homeostasis and energy balance [13,14,15]. Numerous studies have shown that gut microbes assist in coordinating energy balance under cold stress and activate brown adipose tissue thermogenesis [16]. The hindgut serves as a secondary fermentation chamber in ruminants, where approximately 17% of digestible fiber is fermented producing substantial amounts of short-chain fatty acids (SCFAs) [17], including acetate, propionate, and butyrate, significant energy sources for the host [18]. The cecum, being the most anterior section of the hindgut, not only harbors the highest microbial diversity and density in the sheep intestine [19] but is also the primary site for SCFA synthesis in the hindgut [20]. Notably, the cecal epithelial cells express high levels of SCFA transporters (MCT1 and SMCT1), ensuring efficient absorption of these metabolites during periods of energy crisis [21, 22]. SCFAs absorbed in the large intestine account for 8–17% of total absorption in the digestive tract of ruminants, providing 5–10% of metabolic energy [23]. Propionate is a critical precursor for glucose that is synthesized in the liver via gluconeogenesis to support systemic metabolic processes. Furthermore, the cecum contains lymphoid tissues that interact with the microbiota, modulating mucosal immune responses via antigen presentation and helping to suppress inflammatory energy expenditure [24]. While the specific mechanisms of cecal function under cold stress remain unclear, its central role in SCFA metabolism, microbiota-host interactions, immune regulation, and energy balance has been well documented [25,26,27,28]. Therefore, understanding how the cecum regulates cold tolerance through microbial and metabolic pathways offers a promising avenue for uncovering the adaptive strategies of ruminants to cold stress.

Sheep from different growth environments develop distinct adaptive traits [29,30,31]. Hulunbuir sheep and Hu sheep, two ecotypes from different regions of China, have adapted to contrasting climates. Hulunbuir sheep inhabit the Hulunbuir grasslands in eastern Inner Mongolia, a mid-temperate, semiarid continental climate zone with an annual average temperature of – 4.4 °C–3.2 °C and winter temperatures of approximately − 25 °C. They are primarily raised through free grazing and have developed high roughage tolerance and strong cold resistance due to the prolonged period of withered grass. By contrast, Hu sheep, native to the Taihu Basin, live in a subtropical monsoon climate with an average temperature of 15–17 °C and winter temperatures from 5.0 to 12.0 °C. Mainly reared in pens, Hu sheep adapt well to hot and humid conditions. Previous studies indicate that Hu sheep experience a significant decrease in average daily gain (ADG) under cold stress, and Hulunbuir sheep do not, suggesting breed-specific cold adaptation mechanisms [118].

To maintain energy balance during cold exposure, animals typically increase their feed intake to boost energy demands. Herein, rationed feeding was implemented to control for intake and appetite effects, allowing for a detailed examination of adaptive mechanisms in sheep from different genetic backgrounds and selection pressures. By simulating a cold winter environment, we conducted a comprehensive analysis of gut microbiota, SCFAs, cytokines, and transcriptomes to compare the cold adaptation mechanisms of Hulunbuir and Hu sheep.

Methods

Ethics approval statement

The experimental procedures, animal handling, and sample collections were approved by the Animal Welfare and Experimental Ethics Committee of the Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences (protocol number: CAS201810082). In addition, the current study’s experimental procedures were in accordance with ARRIVE guidelines [32].

Experimental design and sample collection

The experimental design and administration were consistent with the previous study [118]. A total of 24 eves (34.65 ± 2.20 kg; mean ± SD) from Hu and Hulunbuir breeds were divided into four groups: Hu sheep at room temperature (HuRT), Hulunbuir sheep at room temperature (HbRT), cold-stressed Hu sheep (HuC), and cold-stressed Hulunbuir sheep (HbC). Sheep in the cold-stress group were first acclimatized at 15 °C for 1 week, then gradually cooled to – 20 °C over 14 days, and maintained at − 20 °C for 38 days. In the room temperature (RT) group, ambient temperature was maintained at approximately 15 °C. Each sheep received a uniform daily ration of 1200 g (Table S1). The study was conducted at the Gaolan Integrated Experimental Station of Ecology and Agriculture (36°13″ N, 103° 47″ E, 1780 m above sea level) [33] in October 2021. Each sheep was housed in individual metabolic cages (1.2 m × 0.6 m × 1.8 m), providing sufficient space and environmental enrichment. The cecal tissues and contents were collected postmortem via slaughtering after 52 days.

Measurement of SCFA and cytokine concentrations

The determination of SCFAs followed previously reported methods [34]. Briefly, cecal samples were thawed, mixed, and weighed before being combined with ultrapure water. After standing, the mixture was centrifuged, and the supernatant was treated with phosphoric acid and crotonic acid solution. The mixture was incubated, centrifuged again, combined with methanol, and filtered for analysis. The gas chromatography parameters were set as follows: injection port temperature at 250 °C, injection volume of 1 μL, split ratio of 5:1, high-purity N₂ as carrier gas, and a capillary column (DB-FFAP, 30 m × 0.32 mm × 0.25 μm, Agilent, USA). The FID detector temperature was 250 °C, with an H₂ flow rate of 30 mL/min, zero-grade air flow rate of 300 mL/min, and high-purity N₂ as the makeup gas at 20 mL/min. The column oven temperature was programmed to increase from 50 to 190 °C at 25 °C/min, held for 3 min, then increased at 30 °C/min to 400 °C, and held for 1 min for chromatographic detection. Cecal levels of secretory IgA (SIgA), IgG, IL-1β, IL-6, IL-10, and TNF-α were quantified using commercial ELISA kits (BYE80034, BYE80027, BYE80020, BYE80010, BYE80213, and BYE80009, Shanghai Bangyi, Co., Ltd., China). Absorbance was measured at 450 nm using a microplate reader (Spectra MAx 190, Molecular Devices, USA).

Metagenomic sequencing and analyses

Total microbial DNA was extracted using the QIAamp PowerFecal Pro DNA Kit (51,804, QIAGEN). The DNA concentration was measured using a Qubit® dsDNA Assay Kit (Q32851, Thermo Fisher Scientific, Carlsbad, CA, USA) and Qubit® 2.0 Fluorometer (Thermo Fisher Scientific, Carlsbad, CA, USA), and the DNA quality was monitored on 1% agarose gels. Sequencing libraries were generated using the NEB Next® Ultra™ DNA Library Prep Kit (E7370L, NEB). DNA samples were fragmented by sonication to 350 bp, which were end-polished, A-tailed, and ligated. The polymerase chain reaction (PCR) amplification was conducted using Universal Illumina adapter primers P5 (5′-AATGATACGGCGACCACCGAG-3′) and P7 (5′-CAAGCAGAAGACGGCATACGAG-3′), with thermal cycling parameters consisting of an initial denaturation phase at 98 °C for 30 s to ensure complete DNA strand separation, followed by 15 amplification cycles comprising three-stage temperature transitions: 10 s of denaturation at 98 °C to maintain template single-stranded, 30 s of annealing at 65 °C optimized for primer-template binding specificity, and 30 s of extension at 72 °C allowing for Taq polymerase-mediated DNA synthesis, culminating in a terminal elongation step at 72 °C for 5 min to ensure complete extension of all PCR products while minimizing partial duplex formation. The clustering of the index-coded samples was performed on a cBot Cluster Generation System and then sequenced on an Illumina Novaseq 6000 platform by Novogene Co., Ltd. (Tianjin, China). Data underwent quality filtering and trimming using Fastp (v0.23.2) [35]. To remove host contamination, Minimap2 (v2.24-r1122) [36] was employed to align the sequences against the sheep genome (ARS-UI_Ramb_v2.0, GCA_016772045.1). Assembly and contig construction were performed using Megahit (v1.1.2) [37], retaining alleles with a minimum length of 300 bp. Gene prediction was conducted using Prodigal (v2.6.3) [38], generating gene and protein sequence files. Transcript expression was quantified using Salmon software [39], which provides fast and bias-aware quantification. Redundancy in protein sequences was reduced using the cluster module of MMseqs2 to obtain a nonredundant protein set. Functional annotation and analysis of gene functions across all samples were achieved by comparing the nonredundant protein sequences with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using the search module of MMseqs2. In addition, microbial taxonomic annotation of contig sequences was performed using the 2bLCA (dual BLAST-based LCA) algorithm, based on MMseqs2.

Cecal RNA sequencing and analyses

In this investigation, total RNA was extracted from animal specimens using TRIzol Reagent (Thermo Fisher Scientific). Primary evaluation of RNA degradation and contamination was performed through 1% agarose gel electrophoresis. Subsequent verification of RNA purity was achieved by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). Quantitative determination of RNA concentration was carried out with the Qubit RNA HS Assay Kit (Q32855, Thermo Fisher Scientific) on a Qubit® 2.0 Fluorometer (Thermo Fisher Scientific, Wilmington, DE, USA). Finally, RNA integrity was systematically evaluated using the RNA Nano 6000 Assay Kit in conjunction with the Agilent 2100 Bioanalyzer System (Agilent Technologies, Santa Clara, CA, USA). The mRNA sequencing library was prepared using the Hieff NGS Ultima Dual-mode mRNA Library Prep Kit (Yeasen Biotechnology, Shanghai, China), encompassing mRNA extraction, cDNA synthesis, end repair, adenylation, adapter ligation, and library purification. Index tags were incorporated through PCR, with the PCR products subsequently purified using AMPure XP beads (Beckman Coulter, Beverly, USA). Library quality was assessed using Agilent Bioanalyzer 2100 before performing 150 bp paired-end sequencing on the Illumina Novaseq 6000 platform. Raw sequencing data were quality controlled using Fastp software (v0.23.4) [35] by trimming read segments, filtering low-quality sequences, removing junctions, etc., with a minimum length threshold of 36 bp and a minimum quality threshold of 20. Subsequently, sequence comparisons were then performed by referring to the NCBI database Sheep Reference Genome ARS-UI_Ramb_v3.0 (GCA_016772045.2) using STAR software (v2.7.11b) [40]. Ultimately, gene expression count reads were extracted from the STAR-generated BAM files using FeatureCounts (v2.0.6) [41] and then normalized to gene length to calculate transcripts per million. A differential expression analysis was conducted for both groups using DESeq2 [42]. Differentially expressed genes (DEGs) were identified based on a P-value < 0.05 and |fold change|> 2. Constrained principal coordinate analysis (CPCoA) used the ImageGP visualization platform [43]. The enrichment of DEGs in the KEGG pathway was analyzed using ClusterProfiler software [44]. The DEGs were aligned with the STRING database (http://string-db.org/, Version: 12.0) using the R package “stringdb” to analyze protein–protein interactions (PPIs). The PPIs of these DEGs were then visualized in Cytoscape [45]. The top 10 hub genes in the PPI network were computationally screened based on the maximum clique centrality (MCC) method using the CytoHubba plugin [46].

Statistical analysis

The differences in SCFAs and inflammatory cytokines between groups were analyzed using GraphPad Prism 9.5. Omics data analysis was performed using R version 4.3.3. Microbial alpha diversity was assessed using the “microbiotaprocess” package [47] based on Wilcoxon test differences. Beta diversity analysis was conducted using the “vegan” package, with principal coordinate analysis (PCoA) calculated based on Bray–Curtis distances and tested using Adonis. Differential microbes (DEMs) and KEGG pathways were analyzed using the linear discriminant analysis effect size (LEfSe) method from the “microeco” package. MetaCyc pathways were assessed using DESeq2 based on counts, and co-occurrence networks were constructed using the SparCC method from the “microeco” package [48]. Module analysis was performed using “cluster fast greedy,” while zipi analysis and network attributes were assessed, with nodes selected at P < 0.05 and |r|> 0.6. Network visualization was conducted using Gephi. Microbial community assembly analysis utilized the “pctax” package (https://github.com/Asa12138/pctax), incorporating Stegen methods (βNTI: β–nearest taxon index, RCbray–based: Bray–Curtis-based Raup–Crick index), phylogenetic normalized stochasticity ratio (pNST), and Solan neutral community model (NCM) methods. The correlation analysis for various results was conducted using the corr.test function from the “psych” package to calculate Spearman correlation coefficients, with Holm correction applied to the P-values. Variables with P < 0.05 and |r|> 0.7 were retained. Network graphs were constructed using the “igraph” R package and visualized with Cytoscape (version 3.10.2). The partial least squares path modeling (PLS-PM) for temperature, ADG, SCFAs, cytokines, microbial composition, microbial KEGG functions, and transcriptomics was built using the “plspm” R package.

Results

Cold exposure enhances SCFA metabolism and cytokine production in Hu and Hulunbuir sheep

Initially, we analyzed SCFAs in the cecum. Results showed the SCFA composition primarily comprising acetate, propionate, and butyrate (Fig. 1A). Compared with RT, cold exposure in Hu sheep resulted in higher concentrations of acetate, butyrate, and total SCFAs (P < 0.05). In Hulunbuir sheep, cold exposure increased propionate and butyrate levels (P < 0.05), with a trend toward increased total SCFAs (0.05 < P < 0.1). In addition, under cold stress, Hu sheep had lower butyrate concentrations than Hulunbuir sheep (P < 0.05), and the acetate/propionate (A/P) ratio showed no differences among groups (P > 0.05) (Fig. 1B). Upon evaluating cytokines in cecal tissue, we found no significant differences in IgG and IL-6 levels among the groups (P > 0.05). Compared with RT conditions, both breeds exhibited increased IL-1β concentrations during cold stress (P < 0.05). In Hulunbuir sheep, SIgA levels also increased (P < 0.05), and Hu sheep showed a trend toward an increase (0.05 < P < 0.1). However, TNF-α concentrations decreased in Hu sheep (P < 0.05), with a trend toward increased IL-10 levels (0.05 < P < 0.1), whereas no significant changes were observed in Hulunbuir sheep. In addition, at RT, IL-10 concentrations were higher in Hulunbuir sheep compared with Hu sheep (P < 0.05), but no differences were observed after cold exposure (Fig. 1C).

Fig. 1

Effects of cold stress on the cecal SCFAs and cytokines of Hu sheep and Hulunbuir sheep. A Composition of SCFAs. B Differential analysis of major cecal SCFA indices; A/P: acetate/propionate ratio. C Differential analysis of cecal inflammatory cytokines. Differences between two groups were identified using independent two-way ANOVA. Significance is indicated by symbols (**, P < 0.01; *, P < 0.05; #, 0.05 ≤ P < 0.1)

Variations in cecal microbiota diversity, composition, and functions

The metagenomic sequencing process yielded a total of 1,059,181,035 reads, with an average of 44,132,543 ± 2,446,331 reads per sample (mean ± SEM). After quality control and host gene elimination, 1,013,999,112 reads were retained, averaging 42,249,963 ± 2,927,436 reads per sample. Subsequently, de novo assembly produced 6,968,420 contigs with an N50 length of 1866 ± 418 bp and an average of 290,351 ± 50,860 contigs per sample (Table S2).

Metagenomic sequencing was employed to analyze cecal microbial diversity and taxonomy across four groups. The Shannon index (P = 0.009) and Chao1 index (P = 0.002) were significantly reduced in Hulunbuir sheep, with the former showing lower diversity than Hu sheep (P = 0.004) at RT (Fig. 2A, B). Bray–Curtis PCoA analysis revealed significant differences attributed to both temperature (R2 = 0.1, P = 0.005) and breed (R2 = 0.08, P = 0.02). The Adonis test indicated significant differences in Hulunbuir sheep between cold stress and normal temperature conditions, as well as between the two breeds under normal temperatures (Fig. 2C; Table S3). Cecal microbial sequencing uncovered alterations in species abundance at the phylum level (Fig. 2D). To compare differences in hindgut microbiota induced by cold exposure, differential species analysis was conducted using the LEfSe method. The results indicated that in Hu sheep, most microbiota that increased under cold stress belonged to the family Lachnospiraceae (10/16), including Faecalicatena contorta, Ruminococcus gnavus, Roseburia hominis, Anaerostipes, Blautia, Coprococcus, and Dorea. In addition, Pseudoflavonifractor sp. MSJ 30 and Bacteroides xylanisolvens exhibited increased abundance (Fig. 2E). Conversely, the abundance of Eukaryota (7/15), primarily Blastocystis sp. subtype 4, and Tenericutes (4/15), including Acholeplasma sp. CAG:878, as well as Clostridium sp. CAG:594, Clostridium sp. CAG:433, and Lentisphaerae bacterium ADurb.Bin242, significantly decreased. In contrast, in Hulunbuir sheep, the abundance of bacteria (34/36), particularly within Spirochaetes (10/36), such as Treponema, Treponema bryantii, and Treponema rectale, significantly increased. Members of Lachnospiraceae (10/36), such as Roseburia sp. 499 and Ruminococcus gnavus, as well as Prevotella copri, Parabacteroides merdae, Clostridium sp. CAG:510, and Methanosphaera sp. SHI1033, also exhibited increased abundance. Conversely, 26 species exhibited decreased abundance, including Clostridiales bacterium, Sarcina sp. DSM 11001, Corallococcus sp. CAG:1435, Anaerotruncus sp. CAG:390, Eubacterium sp. CAG:180, Methanobrevibacter olleyae, and Methanobrevibacter ruminantium (Fig. 2F). Notably, the DEMs in the appendix of Hu and Hulunbuir sheep were less pronounced under cold stimulation, with a higher abundance of bacterium P3, bacterium F082, and Firmicutes bacterium CAG:124 in Hu sheep, and Lachnospiraceae bacterium in Hulunbuir sheep (Fig. S1A). In addition, at RT, Hulunbuir sheep exhibited higher abundances of Treponema bryantii, Treponema rectale, Prevotella sp. ne3005, Prevotella ruminicola, Prevotella copri, Prevotella sp. tc2-28, and others (Fig. S1B).

Fig. 2

Analysis of changes in cecal microbial diversity and composition under cold stress in Hu and Hulunbuir sheep. A, B Alpha diversity was measured by Shannon and Chao1 indices using the Wilcoxon test. C the PCoA plot based on the Bray–Curtis distance and Adonis method. D The top 10 microbiota at the phylum level for each group. Significant changes in the cecal microbiota of Hu (E) and Hulunbuir sheep (F) were identified using LEfSe (LDA > 2, P < 0.05). G Screening of differential KEGG pathways based on the LEfSe approach (LDA > 2, P < 0.05)

Overall, the cecal microbial functions of Hu and Hulunbuir sheep underwent significant changes following cold stimulation, with no notable differences observed between normal and cold environments (Fig. S2; Table S4). Under cold stress, both breed upregulated pathways related to metabolism, environmental information processing, and cellular processes. Hu sheep exhibited an increase in carbohydrate metabolism (starch and sucrose metabolism, galactose metabolism, pentose and glucuronate interconversions, and amino sugar and nucleotide sugar metabolism), amino acid metabolism, metabolism of other amino acids, and methane metabolism, and pathways associated with human diseases, genetic information processing, and organismal systems were downregulated. Hulunbuir sheep upregulated pathways in amino acid metabolism (valine, leucine, and isoleucine biosynthesis; alanine, aspartate, and glutamate metabolism; and phenylalanine, tyrosine, and tryptophan biosynthesis), carbohydrate metabolism (C5-branched dibasic acid metabolism), energy metabolism (oxidative phosphorylation, and nitrogen metabolism), and pantothenate and CoA biosynthesis, and downregulated pathways in genetic information processing, carbon fixation in prokaryotes, and steroid hormone and carotenoid biosynthesis. Comparative analysis under cold stress revealed that Hulunbuir sheep were enriched in C5-branched dibasic acid metabolism, carbohydrate metabolism, antigen processing and presentation, and retinol metabolism (Fig. 2G). Conversely, Hu sheep were enriched in the metabolism of other amino acids, terpenoids and polyketides, fructose and mannose metabolism, and ubiquinone and other terpenoid-quinone biosynthesis (Fig. S1C). At RT, Hulunbuir sheep exhibited enrichment in carbohydrate metabolism, pentose and glucuronate interconversions, phenylalanine metabolism, and methane metabolism, and these pathways were reduced in Hu sheep (Fig. S1D).

We constructed gut microbial networks to capture intermicrobial associations, illustrating the complex relationships and interdependencies within microbial communities. By selecting the top > 0.01% of species in each group and calculating correlations using the SparCC method to construct co-occurrence networks (Fig. 3A), we observed an increase in edges in the networks of Hu sheep and Hulunbuir sheep after cold stimulation (Hu sheep 7769 > 5538; Hulunbuir sheep 7076 > 5947), and the number of nodes remained similar. Modular analysis was performed using the “cluster fast greedy” method, analyzing modular connectivity (Zi) to identify network and module hubs. Results indicated that in the HuC network, there were six module hubs (Alistipes communis, Roseburia sp. 499, Galactobacillus timonensis) and nine network hubs (Bacteroides sp. CAG:545, Unclassified_g_Acetatifactor, and Clostridium sp. CAG:628). In the HuRT network, there were four module hubs (Bacterium P201, Coprobacillus sp. CAG:605, and Prevotella brevis) and one network hub (Bacteroides fragilis). In the HbC network, there were 6 module hubs (Muribaculaceae bacterium Isolate − 013, Clostridium sp. CAG:632, Parabacteroides distasonis) and 6 network hubs (Blastocystis sp. subtype 1, Ruminococcus sp. HUN007, and Clostridium sp. CAG:510). In the HbRT network, there were two module hubs (Ruminococcus sp. CAG:177, and Negativibacillus massiliensis) and six network hubs (Prevotella copri, Unclassified_g_Roseburia, and Firmicutes bacterium CAG:555) (Fig. 3B). After cold stimulation, the microbial co-occurrence network properties of both Hu and Hulunbuir sheep exhibited several similar changes. Compared with the normal-temperature group, Hu sheep exhibited significant changes in average degree, average path length, clustering coefficient, density, heterogeneity, centralization, and modularity. By contrast, these properties changed less in Hulunbuir sheep, but the proportion of negative correlations increased (HbC 48.16%; HbRT 46.01%). In addition, the vulnerability of the microbial network increased in Hu sheep (HuC 0.305; HuRT 0.278), while no change was noted in Hulunbuir sheep (HbC 0.296; HbRT 0.291) (Fig. S3).

Fig. 3

Construction of the co-occurrence network of cecal microbiota using the SparCC method. A Microbial co-occurrence network was constructed for different groups. Node colors represent microbes at different phylum levels, red lines indicate positive correlations, and green lines indicate negative correlations. B Network intramodule connectivity (Zi) and intermodule connectivity (Pi) for calculating network and module hubs

Understanding the roles of deterministic and stochastic processes in the assembly of microbial communities aids in elucidating ecological processes and community responses to environmental changes. We employed a comprehensive approach using Stegen (βNTI and RCbray-based), pNST, and Solan NCM to study the internal driving forces shaping specific community types. Results indicated that at normal temperatures, homogenizing dispersal (Hu sheep 66.67%; Hulunbuir sheep 60%) was more significant than heterogeneous selection (Hu sheep 33.33%; Hulunbuir sheep 40%) in the microbial communities of both sheep breeds. After cold stimulation, homogenizing dispersal in Hu sheep increased (73.33% > 66.67%), and heterogeneous selection decreased (26.67% < 33.33%). Conversely, in Hulunbuir sheep, the contribution of homogenizing dispersal decreased (20% < 40%), and heterogeneous selection increased and became dominant (80% > 40%) (Fig. 4A), indicating that microbial community composition is primarily driven by environmental selection. pNST assesses the relative importance of deterministic and stochastic processes in community construction. The results in Fig. 4B show that the pNST values for each group were > 0.5, indicating that stochastic processes were dominant. The pNST of HbC was significantly lower than that of HbRT and HuC, indicating a reduced contribution of stochastic processes. The Solan NCM results (Fig. 4C) showed that the overall fit of the NCM was similar among the groups, with R2 values for HuC (0.615) and HuRT (0.601) being close, and R² of HbC (0.609) slightly lower than that of HbRT (0.621). After cold stimulation, the construction of cecal microbial communities in Hulunbuir sheep was less influenced by stochastic processes. The Nm of HuC was higher than that of HuRT (1,052,866 > 986,436), and that of HbC was lower than that of HbRT (862,667 < 14,138,104), indicating different diffusion limitations on the cecal microbiota of Hu and Hulunbuir sheep after cold stimulation. Cold stress reduced the influence of neutral processes in the microbial communities of Hulunbuir sheep, imposing greater constraints on microbial dispersal, whereas Hu sheep were less affected.

Fig. 4

Effect of cold stress on microbial community assembly processes in the cecum of Hu and Hulunbuir sheep. A Categorization of deterministic and stochastic processes underlying cecal microbiota diversity using βNTI and RCBray metrics. Hete_S: heterogeneous selection, βNTI < − 2; Homo_S: homogeneous selection, βNTI > 2; D_limit: dispersal limitation, |βNTI|< 2 and RCBray > 0.95; Homo_D: homogenizing dispersal, |βNTI|< 2 and RCBray < − 0.95; Undominated: |βNTI|< 2 and |RCBray|< 0.95 along the sheep cecal microbiota. βNTI, β-nearest taxon index; RCBray, Bray–Curtis-based Raup–Crick Index. B Quantifying the stochasticity in ecological processes using phylogenetic normalized stochasticity ratio (pNST). C Fitting the neutral model to each group of cecal microbiota. R.²: goodness of fit for the overall neutral community model, Nm: product of metacommunity size (N) and migration rate (m) (Nm = N*m)

Cold stress alters gene transcript expression in cecal tissues of Hu and Hulunbuir sheep

Further transcriptome sequencing of cecal tissues revealed that 13,633 genes were expressed. DEGs were identified using thresholds of |Log2FoldChange|> 1 and P < 0.05 (Fig. 5A). In the HuC vs. HuRT comparison, 543 DEGs were identified: 362 upregulated and 181 downregulated DEGs. In the HbC vs. HbRT comparison, 913 DEGs were identified: 361 upregulated and 552 downregulated DEGs. In the HbC vs. HuC comparison, 720 DEGs were identified: 340 upregulated and 380 downregulated DEGs. In the HbRT vs. HuRT comparison, 514 DEGs were identified: 193 upregulated and 321 downregulated DEGs (Fig. 5B, C).

Fig. 5

Effect of cold stress on the cecal transcriptional profile of Hu and Hulunbuir sheep. A CPCoA analysis. B DEG volcano plot. C DEG Venn plot. D KEGG pathway analysis for the DEGs of HuC vs. HuRT. E network diagram of KEGG pathway and DEGs correspondence in HuC vs. HuRT. F KEGG pathways analysis for DEGs of HbC vs. HbRT. G Network diagram of the KEGG pathways and DEGs correspondence in HbC vs. HbRT. H and I Top 20 hub genes were screened from the PPI networks of the DEGs using the MCC algorithm of the CytoHubba plugin. The color intensity indicates the MCC rating, with red indicating a higher rating

KEGG pathway enrichment analysis revealed that in the HuC vs. HuRT comparison, DEGs were enriched in 24 pathways. These included 10 pathways related to organismal systems (e.g., complement and coagulation cascades, thermogenesis, and chemokine signaling), 8 pathways associated with human diseases (e.g., prion disease), 3 pathways related to environmental information processing (e.g., viral protein interaction with cytokines), and pathways involved in energy metabolism (oxidative phosphorylation), and ribosome function (Fig. 5D). Characterization of DEGs in specific KEGG pathways showed that most genes were upregulated (fold change > 2). Notably, genes such as CCL1, CCL5, CCL19, CCR10, CCR7, CX3CR1, CXCL13, LOC100101238, and LOC101119572 were enriched in pathways involving viral protein interaction with cytokine receptors, cytokine–cytokine receptor interaction, and chemokine signaling. Genes such as ATP5ME, NDUFAF7, LOC101121518, LOC101116843, and LOC101119721 were enriched in thermogenesis and oxidative phosphorylation (Fig. 5E). GSEA analysis also indicated upregulation of these pathways and various disease-related pathways (Fig. S4A).

In the HbC vs. HbRT comparison, seven pathways were enriched, including three related to the immune system (e.g., complement and coagulation cascades, PPAR signaling, and IL-17 signaling), two related to signaling molecules and interactions (e.g., viral protein interaction with cytokines), one related to human diseases (COVID-19), and cytoskeleton in muscle cells (Fig. 5F). GSEA analysis also showed upregulation of these pathways and various disease-related pathways (Fig. S4B). IL6 and IL6R were enriched in pathways involving viral protein interaction with cytokine receptors, cytokine–cytokine receptor interaction, and COVID-19. Genes such as CCR2, CSF1R, CXCL14, CXCR2, LOC101108811, LOC101115044, and LOC101117971 were enriched in pathways involving viral protein interaction with cytokine receptors and cytokine–cytokine receptor interaction (Fig. 5G).

DEG enrichment analysis for HbC vs. HuC identified four immune system pathways and three disease-related pathways, as well as transcription-related pathways such as the ribosome, ribosome biogenesis in eukaryotes, and viral protein interaction with cytokine and cytokine receptors. Genes such as CCL5 and LOC101114285 were enriched in pathways such as viral protein interaction with cytokine and cytokine receptors, toll-like receptor signaling, and cytosolic DNA-sensing. CXCL11, LOC100101238, and LOC101114535 were enriched in viral protein interaction with cytokine and cytokine receptors, and toll-like receptor signaling pathways. IRF7 and LOC114113004 were enriched in toll-like receptor signaling and cytosolic DNA-sensing pathways (Fig. S5A, B; Fig. S4C). In the HbRT vs. HuRT comparison, KEGG enrichment analysis identified six pathways: four related to human diseases and the remaining involved natural killer cell-mediated cytotoxicity and ribosome. Seventeen LOC genes were simultaneously enriched in the ribosome and COVID-19 pathways, and six genes were enriched in graft-versus-host disease and natural killer cell-mediated cytotoxicity (Fig. S5C, D; Fig. S4D).

PPI analysis of DEGs using the MCC algorithm in CytoHubba, identified the top 20 hub genes. In the HuC vs. HuRT comparison, all hub genes were upregulated. Key genes included IL1B, seven chemokine receptor family members (e.g., CCR7, CXCL13, CX3CR1, and CCL1), three complement system components (C1QA, C1QB, and C1QC), and others such as CD74, MRPL40, MRPL23, and NHP2 (Fig. 5H). In the HbC vs. HbRT comparison, six downregulated genes were identified, primarily IL6, NPY, and COL2A1. Additionally, there were 12 upregulated genes, including four from the matrix metalloproteinase (MMP) family (MMP9, MMP13, MMP1, MMP12), as well as CXCR2, IL1RN, GCG, GNB3, and SELL (Fig. 5I).

Integrative analysis of multiomics and phenotypic associations

To explore relationships among cytokines, SCFAs, DEMs, and genes in the cecum, we calculated Spearman correlations, retaining those with |r|> 0.7 and P < 0.05. The correlation network for Hu sheep had fewer nodes (58) and edges (156) compared with that of Hulunbuir sheep, which had 75 nodes and 306 edges. The Hu sheep network had a higher proportion of negative correlations (56.41%), whereas the Hulunbuir sheep network had more positive correlations (55.88%). In Hu sheep, butyrate and total SCFAs were significantly positively correlated with Bacteroides xylanisolvens, bacterium 1xD42-87, CCR7, and C1QB, and significantly negatively correlated with Clostridium sp. CAG:1000, Clostridium sp. CAG:594, and cell cycle yeast. SIgA and IL-10 were significantly positively correlated with Firmicutes bacterium CAG:646 and bacterium 1XD42-54, and significantly negatively correlated with Blastocystis sp. subtype 1, antigen processing and presentation, and parathyroid hormone synthesis, secretion, and action. ADG was negatively correlated with bacterium 1XD42-54, amino sugar and nucleotide sugar metabolism, cyanoamino acid metabolism, and CCL1, and positively correlated with Acholeplasma sp. CAG:878 (Fig. 6A). In Hulunbuir sheep, more complex relationships were observed among multiple indicators: SIgA, butyrate, flagellar assembly, and Firmicutes bacterium ADurb.Bin080 showed numerous correlations. SIgA was positively correlated with propionate, butyrate, Prevotella sp. CAG:1031, Treponema sp. C6A8, Treponema sp. JC4, MMP13, MMP9, C5 branched dibasic acid metabolism, and oxidative phosphorylation, among other parameters. Butyrate was significantly correlated with several microbes and functions, being positively correlated with Parabacteroides merdae, Roseburia sp. 499, Ruminococcus gnavus, Treponema sp. C6A8, C5-branched dibasic acid metabolism, flagellar assembly, and two-component system. Flagellar assembly was significantly positively correlated with multiple microbes, including Parabacteroides merdae, Prevotella sp. CAG:1031, Roseburia sp. 499, Ruminococcus gnavus, and Treponema rectale. None of the indicators showed strong correlations with ADG in Hulunbuir sheep (Fig. 6B). We constructed PLS-PM models for Hu and Hulunbuir sheep to comprehensively analyze potential relationships among phenotypes, transcriptomes, and microbial composition and function. The two breeds exhibited distinct differences. In Hu sheep, complex potential relationships among various indicators were observed, mostly negative but not significant. Only microbial diversity and composition showed a significant negative correlation, and microbial function had a significant positive correlation with SCFAs. Temperature had a somewhat negative impact on cecal gene transcription, and microbial function negatively affected cytokines (Fig. 6C). For Hulunbuir sheep, fewer interactions were observed, resulting in a model with a higher goodness of fit. Temperature had a significant negative impact on microbial diversity and the transcriptome. Microbial diversity influenced microbial composition, which in turn had a significant positive effect on microbial function, SCFAs, and cytokines. In addition, we analyzed the standardized effects of different indicators on ADG and cytokines. In Hu sheep, microbes had a substantial negative effect on ADG, and microbial function had a positive effect. Microbes had a large negative effect on cytokines. In Hulunbuir sheep, SCFAs had a significant negative effect on ADG, and cytokines were positively influenced by microbes and their functions and negatively by temperature and microbial diversity (Fig. 6D).

Fig. 6

Correlation network among cecal SCFAs, cytokines, microbial functions, DEMs, and DEGs in Hu sheep (A) and Hulunbuir sheep (B). The PLS-PM model provided an integrated analysis of the relationships among various phenotypes, transcriptomes, and microbial community functions in both Hu sheep (C) and Hulunbuir sheep (D)

Discussion

Cold exposure affects nutrient digestion, metabolism, and gut microbiota in sheep, negatively affecting their health and growth. Genetic background and growth environment influence the gut microbiota and metabolic preferences of animals. Numerous studies have identified specific microbes and metabolites in different phenotypes that enhance animal performance and health [49,50,51,52,53,54]. Previous experiments found that under cold stress, the ADG of Hu sheep was significantly lower than that of the control group, and Hulunbuir sheep showed no significant change. Herein, we designed a controlled feeding and temperature experiment to investigate the adaptation mechanisms of different sheep breeds to cold exposure. We combined analyses of gut SCFA metabolism, microbial metagenomics, cecal transcriptomics, and cytokine profiling. The results revealed distinct microbial responses, SCFA metabolism, and transcriptional regulation mechanisms between Hu and Hulunbuir sheep. Hulunbuir sheep demonstrated greater microbial plasticity and metabolic shifts, enhancing energy utilization and immune responses to maintain ADG. By contrast, Hu sheep exhibited smaller microbial changes, with adaptations in acetate metabolism and thermogenesis leading to a decrease in ADG.

SCFAs are the most abundant fermentation metabolites in the gastrointestinal tract of ruminants, with acetate, propionate, and butyrate comprising over 95% of total SCFAs [55]. Consistent with the findings of this study, previous research has demonstrated that ovine subjects exhibit an increase in SCFA concentrations within the gastrointestinal tract during colder seasons [56,57,58]. In this study, under restricted caloric feeding, cold exposure promoted SCFA metabolism, leading to increased concentrations of SCFAs, particularly butyrate. Butyrate serves as a primary energy source for intestinal cells and offers multiple benefits, including anti-inflammatory properties, enhanced gut barrier function, and improved mucosal and systemic immunity, thereby reducing infections by gut bacteria or pathogens [59]. In addition, butyrate activates thermogenesis, regulating body temperature and energy expenditure by acting as a signaling molecule and activating GPR43, GPR41, and GPR109a receptors in brown and white adipose tissues [60, 61]. Studies in rats have also found that butyrate can mitigate the adverse effects of cold exposure [62]. In Hu sheep, acetate concentrations increase, and in Hulunbuir sheep, propionate levels increase. Acetate, a primary product of structural carbohydrate (e.g., cellulose) metabolism, predominantly enters the circulatory system and is delivered to tissues such as adipose and mammary glands, where it is converted to acetyl-CoA, supporting the metabolism of fats, proteins, and carbohydrates. Propionate, primarily derived from nonstructural carbohydrate metabolism, absorbs some H2 energy during its production, making propionate fermentation the most energy-efficient among SCFAs [63]. Propionate serves as a crucial energy substrate for ruminants, with glucose required for ruminant metabolism mainly produced through the hepatic gluconeogenesis of propionate. Studies have shown increased hepatic propionate uptake in sheep during cold exposure to synthesize glucose [64]. The elevated cecal propionate in Hulunbuir sheep may reduce energy loss during fermentation and serve as an energy substrate, enhancing their cold adaptation. In mice, increases in both propionate and butyrate have been associated with improved cold adaptation [13]. Compared with Hu sheep, the enhanced propionate and butyrate metabolism in Hulunbuir sheep likely improves environmental adaptability.

Cold exposure also triggers immune responses, altering various inflammatory factors. Exposure to stressors induces the release of pro-inflammatory cytokine IL-1β [65, 66]. Our findings align with this, because cold exposure also led to IL-1β release, intensifying inflammation [67, 68]. SIgA in the sheep gut forms a protective layer on the mucosal surface, preventing pathogen invasion and playing a crucial role in mucosal protection [69, 70]. Cold stress has been shown to increase IgA expression in the intestines of chickens [71], enhancing gut immune function to resist cold environments [72]. IL-10, a cytokine involved in inflammation and immune suppression, activates macrophages to inhibit the expression of inflammatory cytokines such as TNF-α, IL-6, and IL-1, thereby reducing inflammation [73]. Similarly, cold exposure leads to increased IL-10 levels in mice and chickens to combat inflammation [74, 75]. After cold exposure, both Hu and Hulunbuir sheep exhibit inflammatory responses and mitigate inflammation through different mechanisms.

Under cold exposure, the gut microbiota diversity in animals undergoes distinct changes [13, 76, 77]. In Hulunbuir sheep, both alpha and beta diversity are altered, indicating significant microbial community adjustments. Lachnospiraceae, a dominant group in the sheep hindgut, is linked to SCFA metabolism [78]. In the ceca of Hu sheep, several Lachnospiraceae species increase, consistent with studies on cold-stressed rodents [79]. Roseburia hominis, Coprococcus, Blautia, and Pseudoflavonifractor can ferment carbohydrates to produce butyrate [80, 81]. Roseburia hominis is considered a potential probiotic, with its butyrate production helping to maintain intestinal epithelial integrity and immune balance [82, 83]. Faecalicatena contorta and Bacteroides xylanisolvens produce acetate as a fermentation end product [84, 85], enhancing acetate and butyrate metabolism in Hu sheep cecum. There is a notable reduction in Blastocystis, a common eukaryotic gut symbiont in humans and animals, often considered potentially pathogenic and associated with gut infections and inflammation [86, 87]. Our findings show that Blastocystis sp. subtype 1 is significantly negatively correlated with the anti-inflammatory factor IL-10. In Hulunbuir sheep, an increase in various species of Treponema is observed. Treponema, a major gastrointestinal species, ferments to produce acetate and succinate [88, 89], with succinate being a key intermediate in propionate metabolism. Prevotella copri ferments cellulose and participates in the succinate pathway, which is related to propionate production. Ruminococcus gnavus possesses a coenzyme A-dependent propionaldehyde dehydrogenase gene in the propanediol pathway, generating propionate [90]. In addition, there is an increase in Lachnospiraceae species involved in butyrate metabolism. Interestingly, cold exposure also alters archaeal populations in Hulunbuir sheep. Methanobrevibacter ruminantium, Methanobrevibacter millerae, and Methanobrevibacter olleyae, key hydrogenotrophic methanogens [91, 92], decrease in abundance, possibly due to enhanced propionate metabolism capturing H₂ [93]. These microbial changes enhance propionate and butyrate metabolism in Hulunbuir sheep ceca. Functional annotation results indicate that Hu sheep microbiota enhances various carbohydrate metabolisms, promoting SCFA production. Enhanced methane metabolism may be linked to increased acetate concentrations, with glucose metabolism producing acetate and H₂ for methane synthesis [94]. Hulunbuir sheep showed increased amino acid metabolism, with alanine, aspartate, and glutamate metabolized into propionate and butyrate [95], and branched-chain amino acids such as valine, leucine, and isoleucine also contribute to butyrate production [96]. Hulunbuir sheep demonstrate stronger carbohydrate metabolism for SCFA production in both normal and cold environments. Enriched pathways in environmental information processing and cellular processes in both sheep breeds facilitate microbial survival and colonization.

Microbial network construction can simulate interspecies relationships and can assess microbiome stability [97]. In the cecal microbiota co-occurrence networks of both Hu and Hulunbuir sheep under cold exposure, there are more connections, indicating tighter microbial interactions. Increased network complexity may relate to metabolic functions [98]. The Hu sheep cold-stress network shows significant changes, with increased vulnerability, and the Hulunbuir sheep network remains more stable, benefiting host health [99]. Understanding the roles of deterministic and stochastic processes in the assembly of microbial communities aids in elucidating ecological processes and the response of communities to environmental changes [100]. In Hu sheep, the microbiota remains relatively stable under cold stress, with stochastic processes dominating and no significant selective pressures observed, aligning with diversity findings. By contrast, Hulunbuir sheep exhibit significant changes in microbial community structure under cold stress, with increased heterogeneous selection indicating strong environmental selection pressures. The cecal microbiota of Hulunbuir sheep demonstrates greater adaptability [101].

Cold exposure induces transcriptome reprogramming in multiple tissues of animals [102], with distinct changes observed in the cecal transcriptomes of Hu and Hulunbuir sheep. In Hu sheep, cold exposure results in the upregulation of DEGs in the cecum, enhancing various functions. The upregulation of pathways such as oxidative phosphorylation, thermogenesis, and ribosome biogenesis indicates increased energy metabolism and protein synthesis, aiding environmental adaptation [76, 103]. ATP5ME, a subunit of ATP synthase, and genes such as LOC101117181, LOC101119721, LOC101121036, LOC101121518, LOC114110252, LOC114115573, and LOC121817235 are involved in the assembly of cytochrome c oxidase. NDUFA7 and NDUFB8 are components of mitochondrial electron transport chain complex I. The upregulation of these genes promotes ATP synthesis, enhancing the energy production necessary for cold adaptation. Complement and coagulation cascades, chemokine signaling pathways, viral protein interactions with cytokines and cytokine receptors, and cytokine-cytokine receptor interactions are upregulated in response to cold exposure in both the Hu sheep and Hulunbuir sheep groups. This suggests that cold exposure may enhance immune defense and mitigate inflammation through cytokine regulation [104, 105]. The upregulation of IL1B is identified as a key gene in Hu sheep under cold stress, aligning with cytokine assay results. CCR7, CXCL13, CX3CR1, and CCL1 belong to the chemokine receptor family. In Hulunbuir sheep, the core genes under cold stress are mainly from the MMP family, including MMP9, MMP13, MMP1, and MMP12, which facilitate tissue repair and regulate proinflammatory factors TNF-α and IL-6, thereby enhancing immune defense mechanisms and improving pathogen clearance [106]. Compared with Hu sheep, Hulunbuir sheep exhibit downregulation in thermogenesis, oxidative phosphorylation, and various disease-related pathways under cold stress. Hulunbuir sheep do not activate multiple energy metabolism and inflammatory pathways and reduce protein synthesis, which may be related to enhanced propionate metabolism to conserve energy. Cecal tissues show increased transcription for inflammation resistance.

The phenotype and multiomics correlation network elucidate the interrelationships and impact magnitudes of various indicators. In both Hu and Hulunbuir sheep networks, butyrate occupies a central position, linking specific microbes with functions and tissue transcription, consistent with recent studies highlighting butyrate’s role in cold adaptation [60, 62, 107]. Hulunbuir sheep exhibit significant adjustments in microbial diversity and community assembly, with complex relationships among various indicators, similar to the PLS-PM results. Hulunbuir sheep significantly alter cecal microbial diversity, composition, function, and SCFA metabolites to enhance energy utilization and cold adaptation, maintaining similar ADG changes compared with normal temperatures.

Environmental stress poses a major challenge for sustainable livestock development, and enhancing cold adaptability in sheep holds significant potential for healthy growth. This study revealed breed-specific differences in cecal microbiota, SCFA metabolism, and cecal transcriptional responses under cold stress. Distinct genetic backgrounds shaped divergent microbial metabolic profiles and host tissue adaptations. Notably, cecal propionate and butyrate, along with their associated microbial consortia, may play pivotal roles in cold adaptation. While previous studies have employed environmental simulations (controlled temperature, wind speed, etc.) [4, 108,109,110] housing environments [111,112,113], or seasonal comparisons [114,115,116,117], these approaches often introduce confounding variables. This study implemented strict temperature control and restricted feeding regimens in two short-tailed ewe breeds, and future investigations can expand the cohort size while systematically investigating cold adaptation mechanisms in rams, wethers, and periparturient ewes. Furthermore, although we observed different responses to specific microbes, propionate, and butyrate, which have been extensively characterized in mouse models [14,15,16,17], validation in sheep remains essential to elucidate the interactions between gut microbiota, SCFAs, and host genetics for energy metabolism. We will progressively validate these findings through in vitro ovine cell models followed by in vivo experiments.

Conclusions

In summary, Hu sheep exhibit an increased abundance of Roseburia hominis, Coprococcus, and Faecalicatena contorta in the cecum, enhancing acetate and butyrate fermentation, which requires tissue thermogenesis and encounters greater immune challenges. By contrast, Hulunbuir sheep display a higher abundance of Prevotella copri, Ruminococcus gnavus, and Roseburia sp. 499, with a reduction in methanogens such as Methanobrevibacter ruminantium, Methanobrevibacter millerae, and Methanobrevibacter olleyae. This shift promotes propionate metabolism for improved energy efficiency and increased butyrate production to strengthen immune response, thereby supporting cold adaptation. Our findings suggest that cecal microbial propionate and butyrate metabolism may contribute to cold adaptation in sheep, with breed-specific strategies observed. However, further validation is required to establish causal relationships between microbial shifts and host physiological outcomes.