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.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s40168-025-02096-9.

   Keywords: Cold stress, Cecal microbiome, Short-chain fatty acids,
   Transcriptome, Immune response

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 [[42]1–[43]6]. In addition, cold stress can lead to various
   diseases, such as respiratory infections [[44]7, [45]8] and intestinal
   inflammations [[46]9–[47]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 [[48]13–[49]15]. Numerous
   studies have shown that gut microbes assist in coordinating energy
   balance under cold stress and activate brown adipose tissue
   thermogenesis [[50]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) [[51]17], including acetate, propionate, and butyrate,
   significant energy sources for the host [[52]18]. The cecum, being the
   most anterior section of the hindgut, not only harbors the highest
   microbial diversity and density in the sheep intestine [[53]19] but is
   also the primary site for SCFA synthesis in the hindgut [[54]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 [[55]21, [56]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 [[57]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 [[58]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 [[59]25–[60]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 [[61]29–[62]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 [[63]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 [[64]32].

Experimental design and sample collection

   The experimental design and administration were consistent with the
   previous study [[65]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)
   [[66]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
   [[67]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[2] 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[2] flow rate of 30 mL/min, zero-grade air flow
   rate of 300 mL/min, and high-purity N[2] 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 ([68]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)
   [[69]35]. To remove host contamination, Minimap2 (v2.24-r1122) [[70]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) [[71]37], retaining alleles with
   a minimum length of 300 bp. Gene prediction was conducted using
   Prodigal (v2.6.3) [[72]38], generating gene and protein sequence files.
   Transcript expression was quantified using Salmon software [[73]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 ([74]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)
   [[75]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) [[76]40]. Ultimately, gene expression count reads were
   extracted from the STAR-generated BAM files using FeatureCounts
   (v2.0.6) [[77]41] and then normalized to gene length to calculate
   transcripts per million. A differential expression analysis was
   conducted for both groups using DESeq2 [[78]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 [[79]43]. The enrichment of
   DEGs in the KEGG pathway was analyzed using ClusterProfiler software
   [[80]44]. The DEGs were aligned with the STRING database
   ([81]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 [[82]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
   [[83]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 [[84]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 [[85]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
   ([86]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. [87]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. [88]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. [89]1C).

Fig. 1.

   [90]Fig. 1
   [91]Open in a new tab

   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. [92]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. [93]2C; Table S3). Cecal microbial sequencing
   uncovered alterations in species abundance at the phylum level
   (Fig. [94]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. [95]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. [96]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.

   [97]Fig. 2
   [98]Open in a new tab

   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. [99]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. [100]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. [101]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.

   [102]Fig. 3
   [103]Open in a new tab

   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. [104]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. [105]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. [106]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^2 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.

   [107]Fig. 4
   [108]Open in a new tab

   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.^2: 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. [109]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. [110]5B,
   C).

Fig. 5.

   [111]Fig. 5
   [112]Open in a new tab

   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. [113]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. [114]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. [115]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. [116]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. [117]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. [118]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. [119]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. [120]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. [121]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. [122]6D).

Fig. 6.

   [123]Fig. 6
   [124]Open in a new tab

   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