Multi-omics insights into microbiome-rumen epithelium interaction mechanisms underlying subacute rumen acidosis tolerance in dairy goats

To investigate the impact of different proportions of rumen degradable starch on rumen fermentation, continuous pH monitoring was conducted in dairy goats for 14 h after morning feeding. The pH of high-grain whole corn SARA-susceptible group (HGW-SARA) goats remained below 5.8 for over 6 h on both days 56 and 63 of the trial period, surpassing the critical threshold for diagnosing SARA occurrence (P < 0.001, Fig. 1A; Additional file 1: Fig S1A-1B). In contrast, none of the healthy goats, including those in the low grain whole corn group (LGW-CON), high grain whole corn SARA-tolerant group (HGW-Health), and high grain crushed corn SARA-tolerant group (HGC-Health), met the criteria for SARA occurrence (Fig. 1A). Significantly higher concentrations of ruminal isobutyrate and isovalerate were detected in HGW-SARA goats than in LGW-CON goats and HGC-Health goats (P < 0.05, Fig. 1B). Compared with HGW-Health goats, only the concentration of isobutyrate was higher in HGW-SARA goats (P < 0.05, Fig. 1B). There were no significant differences in rumen VFAs ratio, acetate to propionate ratio, and other VFA concentrations among the four groups (P > 0.05, Fig. 1B; Additional file 1: Fig S2A-S2B). Additionally, when comparing the differences in NH₃-N, and lipopolysaccharide (LPS) in the rumen and LPS in plasma, results showed that NH₃-N concentration was significantly higher in HGW-SARA goats comparing to LGW-CON goats (P < 0.05), while there were no remarkable differences observed in the concentration of LPS within either the rumen or plasma samples (P > 0.05, Fig. 1C).

The HE-stained section of the rumen epithelial papilla and the statistics on epithelial thickness revealed significant shedding of stratum corneum cells in HGW-SARA goats (Fig. 1D). In comparison to LGW-CON goats, there were significant decreases in the thickness of stratum granulosum and total epithelium in HGW-Health and HGC-Health goats (P < 0.05) with tighter cell arrangement (Fig. 1D). Following the consumption of high-RDS diets, there was a reduction in the thicknesses of the stratum granulosum, stratum spinosum/basale, and total epithelium to varying degrees (P < 0.05, Fig. 1E).

The rumen microbiome metagenomic data revealed that the overall bacteria in rumen PcoA showed a significant separation between the cluster of HGW-SARA goats and LGW-CON goats (R = 0.24, P = 0.047). However, there was no significant difference in β diversity between HGW-SARA and high grain feeding healthy goats (Fig. 2A). LEfSe analysis revealed enrichment of 16 bacterial species in LGW-CON goats, notably belonging to genus Prevotella and phylum Bacteroidetes (LDA > 3, P < 0.05). On the other hand, 12 bacterial species were enriched in HGW-SARA goats with eight of which belonged to the phylum Proteobacteria (LDA > 3, P < 0.05), as shown in Additional file 1: Fig S3A. In comparison to LGW-CON goats, the relative abundance of GH24, GH25, and GH108 was significantly increased in HGW-SARA goats (P < 0.05, Additional file 1: Fig S3B). Ruminococcus sp, Ruminococcus bromii, and Clostridiaceae bacterium were highly correlated with GH24 and GH25 in HGW-SARA goats (Additional file 2: Table S1; Additional file 1: Fig S3C).

LEfSe analysis was conducted to compare differential bacteria of HGW-SARA and HGW-Health goats, and 31 differential bacteria was screened in LDA > 2, P < 0.05. Verrucomicrobia bacterium, Selenomonadaceae bacterium DSM 108025, Mitsuokella sp AF21-1AC, and Bifidobacterium adolescentis were enriched in HGW-Health goats. Comparing with HGW-SARA goats, unclassified-g-Parabacteroides was enriched in HGC-Health goats (Fig. 2B). It was found that only Proteobacteria and Fusobacteria of differential bacteria selected in species level were not enriched in the rumen of healthy dairy goats with high RDS (Additional file 1: Fig S4A). Comparing to HGW-SARA goats, GH43_31, GH30 and GH13_30 were significantly enriched in HGW-Health goats (LDA > 2, P < 0.05), while GH30_5 (LDA = 2.01, P = 0.028) was significantly enriched in HGC-Health goats (Fig. 2C). Phenylalanine, tyrosine and tryptophan biosynthesis, valine, leucine and isoleucine biosynthesis, C5-Branched dibasic acid metabolism, and others were significantly enriched in HGW-Health goats compared with HGW-SARA goats (LDA > 2, P < 0.05, Fig. 2D). However, there were no remarkable differences in the level 3 KEGG pathways between HGC-Health and HGW-SARA goats.

To further delve into the function and significance of rumen microbiome enriched in HGW-Health, we conducted a comparison of the variances in KEGG enzymes between HGW-Health and HGW-SARA goats. Additionally, a total of 56 KEGG enzymes were found to be enriched (LDA > 2, P < 0.05), with an impressive count of 50 being enriched in HGW-Health goats (Additional file 1: Fig S4B), contributing to essential roles across 29 KEGG pathways. Notably, specific enzymes such as 4.2.1.51 and 2.6.1.57 are involved in phenylalanine, tyrosine and tryptophan biosynthesis, and 4.3.1.19, 4.2.1.35, and 4.2.1.33 participated in valine, leucine, and isoleucine biosynthesis, respectively; and 4.2.1.35 is involved in C5-branched dibasic acid metabolism (Additional file 1: Fig S4C). Correlation analysis revealed significant associations between 18 species of bacteria and differential KEGG pathways, as well as related differential KEGG enzymes. Specifically, Selenomonadaceae bacterium DSM 108025 showed positive correlations with 4.3.1.19, 4.2.1.51, and 2.6.1.57 (rho > 0.6, P < 0.05). Additionally, Mitsuokella sp AF21-1AC was positively associated with the enzyme activity of 2.6.1.57 (rho = 0.867, P = 0.003) and C5-Branched dibasic acid metabolism (rho = 0.644, P = 0.044, Fig. 2E). These findings were summarized in Fig. 2F, which shows that these enzymes primarily participated in amino acid biosynthesis such as tyrosine, L-leucine, L-isoleucine, and L-valine while also regulating Threonine metabolism.

Receiver operating characteristic (ROC) analysis was employed to ascertain whether the 4 species of bacteria enriched in HGW-Health goats could function as biomarkers for distinguishing between SARA-susceptible and SARA-tolerant goats on the same diets. It was revealed that 3 out of these 4 species were able to effectively diagnose and classify SARA-susceptible and SARA-tolerant goats with an AUC ≥ 0.88 (P < 0.05). Specifically, Bifidobacterium adolescentis and Mitsuokella sp AF21-1AC exhibited an AUC = 0.88, while Verrucomicrobia bacterium demonstrated an impressive AUC = 0.92 (Fig. 2G).

To further investigate whether alterations in the rumen microbiota influence the host healthy/unhealthy state, we conducted a comparative analysis of gene expression in the host rumen epithelium. Comparing to HGW-SARA goats, 124 genes exhibited downregulated expression and 159 genes showed upregulated expression in HGW-Health goats (P < 0.05, |log2FC|> 1, Additional file 1: Fig S5A). Meanwhile, when comparing HGC-Health goats to HGW-SARA goats, 123 genes were upregulated and 69 genes were downregulated (P < 0.05, |log2FC|> 1, Additional file 1: Fig S5B). Thirty divergent pathways between HGW-Health and HGW-SARA, as well as 16 differential pathways between HGC-Health and HGW-SARA were enriched (Fig. 3A and B). Notably, 8 pathways were shared across both comparisons (Fig. 3C) including the IL-17 signalling pathway, TNF signalling pathway, PI3K-Akt signalling pathway, and others (P_{HGW-SARA vs HGW-Health} < 0.01, P_{HGW-SARA vs HGC-Health} < 0.05, Additional file 2: Table S2; Fig. 3E). Delving further into the overlap between the two differential comparison sets, 45 genes were consistently differentially expressed (P < 0.05, |log2FC|> 1, Fig. 3D). Of particular interest, genes such as TNFAIP3, MMP1, PTGS2, IL6, and CCL20, which are implicated in the IL-17 signalling pathway, were concurrently downregulated in both HGW-Health and HGC-Health goats relative to HGW-SARA goats (Fig. 3E), highlighting a shared response potentially linked to inflammation and immune regulation.

Single-nucleus RNA sequencing was employed to elucidate the cellular constitution of the rumen epithelium. Thirty three distinct clusters were identified, with cluster 6 and cluster 31 being verified as T helper cell 17 (Th17) populations, and cluster 12 confirmed as Macrophages (Fig. 3F). A comparative assessment of the cellular composition in four groups uncovered an elevated proportion of Th17 cells specifically in HGW-SARA goats (7.17%), a marked increase over LGW-CON (3.29%), HGW-Health (5.25%), and HGW-Health (4.74%). Conversely, the macrophage population displayed a relatively consistent ratio across all four goat types (Fig. 3G). Given that IL-17A and CD163 serve as established biomarkers for Th17 cells and macrophages, respectively, immunofluorescence assays were carried out to map their distribution and quantified their expression levels. The findings highlighted a substantial increase in IL-17A expression in HGW-SARA goats compared to LGW-CON, HGW-Health, and HGC-Health goats (P < 0.001, Fig. 3H). Despite the macrophage ratio remaining uniform across groups, CD163 expression was notably augmented in HGW-SARA goats relative to LGW-CON, HGW-Health, and HGC-Health (P < 0.001), with HGC-Health also exhibiting significantly higher CD163 expression compared to both HGW-Health and LGW-CON (P < 0.001, Fig. 3I).

An in-depth single-nucleus RNA analysis not only revealed the presence of Th17 and Macrophages but also delineated a broad cellular landscape, encompassing various cell types: Epithelial cells across twelve clusters (0, 1, 2, 3, 4, 5, 11, 14, 15, 18, 28, 32), Stem cells in four clusters (7, 8, 17, 21), Fibroblasts in two clusters (9, 16), Smooth muscle cells across five clusters (10, 19, 22, 24, 27), Endothelial cells in four clusters (13, 20, 25, 29), Mesenchymal cells in one cluster (23), and Neurons in two clusters (26, 30) (Fig. 3F). Remarkably, Epithelial cells constituted the majority, with their proportion notably augmented in LGW-CON (53.14%), HGW-Health (50.67%), and HGC-Health (58.10%) compared to HGW-SARA (49.41%) (Fig. 3G). Additionally, a stark contrast emerged in the total cellular count within rumen epithelium. Of note, the Epithelial cell counts nearly doubled in HGW-Health and HGC-Health goats relative to LGW-CON and HGW-SARA (Additional file 2: Table S3). To validate these findings, immunofluorescence staining employing KRT6A, a prominent epithelial marker, disclosed a significant upregulation in HGC-Health goats vis-a-vis HGW-SARA (P < 0.001, Fig. 4A). Delving into the underlying cause of this epithelial cell discrepancy, a subcluster classification of epithelial cells was conducted. The UMAP analysis segregated seven subclusters, with their top three marker genes illustrated in Additional file 1: Fig S6, categorizing epithelial cells into Granule cells (clusters 0, 2, 5, 6), Spinous cells (3, 4), and Basal cells (1) (Fig. 4B). A conspicuous variation in Spinous and Basal cell counts across the groups was discernible. Quantification and proportion analysis confirmed Granule cells as the predominant type, but Basal and Spinous cells exhibited the most substantial numeric disparity, corroborating the UMAP visualization (Fig. 4C; Additional file 2: Table S4).

The expression levels of genes in solute carrier family between HGW-SARA and HGW-Health goats were analyzed, as well as between HGW-SARA and HGC-Health goats. The result showed that the expression levels of certain solute carrier gene families were markedly elevated in the HGW-Health and HGC-Health dairy goats (log₂FC ≥ 1, P < 0.05, Fig. 4D). Reflecting on epithelial barrier function, CLDN4 mRNA expression was significantly enhanced in HGW-Health goats compared to LGW-CON and HGW-SARA (Fig. 4E). Similarly, MCT1, indicative of epithelial absorption capacity, was significantly upregulated in HGW-Health goats (Fig. 4F). Lastly, plasma VFAs assessment revealed that propionate, isobutyrate, and total VFAs concentrations were significantly reduced in HGW-SARA goats compared to LGW-CON and HGW-Health, with no significant variance noted in other VFAs or their ratios among all groups (Fig. 4G; Additional file 1: Fig S7A-B). HGC-Health goats did not exhibit significant deviations in either VFAs concentrations or ratios from LGW-CON and HGW-SARA goats (Additional file 1: Fig S7C-D), underscoring the complex interplay of factors influencing epithelial dynamics and metabolic profiles.

A nontargeted metabolomics approach was adopted to contrast HGW-SARA against its healthy counterparts. A distinct separation between LGW-CON goats and those fed high-RDS diets, albeit without clear differences among HGW-SARA, HGW-Health, and HGC-Health goats (Additional file 1: Fig S8A). A permutation test with a Q²Y intercept less than 0.05 and an R²Y intercept less than 0.3 confirmed the absence of model overfitting, attesting to the model’s robustness and reliability (Additional file 1: Fig S8B).

Compared to those in LGW-CON goats, 37 metabolites were significantly altered (VIP > 1, FDR < 0.05), with 17 elevated and 20 diminished in HGW-SARA goats (Additional file 1: Fig S8C). These alterations were implicated in ten KEGG pathways, separately being alpha-Linolenic acid metabolism, cAMP signalling pathway, renin secretion, regulation of lipolysis in adipocytes, and others (Additional file 2: Table S5; Additional file 1: Fig S8D).

Further comparisons revealed 138 differentially abundant metabolites between HGW-SARA and HGW-Health goats, 107 of which were upregulated in HGW-SARA goats and 31 in HGW-Health goats (VIP > 1, P < 0.05, Fig. 5A). Between the HGW-SARA and HGC-Health groups, 52 differentially abundant metabolites emerged, 39 of which were elevated in the HGW-SARA group and 13 in the HGC-Health group (VIP > 1, P < 0.05, Fig. 5B). Cysteine and methionine metabolism, alpha-Linolenic acid metabolism, Arginine and proline metabolism, Purine metabolism, and Bile secretion were enriched in the comparison between HGW-SARA and HGW-Health (P < 0.05, Fig. 5C). For HGW-SARA versus HGC-Health, alpha-Linolenic acid metabolism, Glutathione metabolism, Pyruvate metabolism, and Thyroid hormone synthesis were significantly highlighted (P < 0.05, Fig. 5D). To delve deeper, KEGG topological analysis illuminated metabolic discrepancies. Between HGW-SARA and HGW-Health, Cysteine and methionine metabolism, valine, leucine, and isoleucine biosynthesis, Arginine and Proline metabolism, and other 12 pathways were distinguished (Fig. 5E). Meanwhile, Pentose phosphate pathway, Pyruvate metabolism, and other 4 pathways holding notable impact in the comparison between HGW-SARA and HGC-Health goats (Fig. 5F).

Correlation analysis was employed to elucidate relationships between these metabolites and key rumen parameters such as fermentation products, the concentration of LPS, NH₃-N, and lactate. The results showed that metabolites elevated in HGW-SARA goats compared to HGW-Health were positively correlated with increased VFAs and NH₃-N concentrations. Notably, 3-methyl pyruvic acid (rho = 0.767, P = 0.001) were significantly positively correlated with the ruminal pH. When HGW-SARA was compared with HGC-Health, the upregulated metabolites in HGW-SARA demonstrated stronger positive associations with lactate and pH and inverse relationships with VFAs levels (Fig. 5G).

To investigate how microbiome and metabolites influence the SARA tolerance of HGW-Health goats, correlations were examined between HGW-Health-enriched microbiome and metabolites integral to metabolic regulation. Bifidobacterium adolescentis emerged as a pivotal microbe, displaying strong positive correlations with metabolites like 3-methyl pyruvic acid (rho = 0.733, P = 0.021), L-Methionine S-oxide (rho = 0.758, P = 0.016), uric acid (rho = 0.745, P = 0.018), and Xanthine (rho > 0.6, P < 0.05, Additional file 2: Table S6; Fig. 5H).

Concomitantly, both metagenomic and metabolomic data pointed to enriched valine, leucine, and isoleucine biosynthesis pathways in the HGW-SARA versus HGW-Health comparison. With a focus on the pivotal metabolite 3-methyl pyruvic acid involved in this pathway and the differential bacterium B. adolescentis, a structural equation model (SEM) was constructed to explore the potential interplay among five key elements: B. adolescentis, 3-methyl pyruvic acid, SLC family genes, MCT1, and epithelial cell proliferation. The model suggested a putative pathway in which B. adolescentis may influence MCT1 expression and mediate the effect of 3-methyl pyruvic acid on epithelial cell proliferation, potentially contributing to the maintenance of epithelial integrity (Fig. 5I).

To negate host-specific influences, an interventional study involving rumen content transplantation was executed to elucidate the roles of microbiome and metabolites in rumen function. Rumen fluids from three HGW-SARA goats were transplanted into HGC-Health goats, labelled S + H, while fluids from three HGC-Health goats were transferred to HGW-SARA goats, designated H + S (Additional file 1: Fig S1B). Monitoring rumen pH at the 7th and 14th day post-transplantation revealed that only the HGW-SARA goats maintained a pH below 5.8 for more than 3 h; in contrast, both the S + H and H + S groups exhibited higher pH values (Fig. 6A). In VFA profiles, S + H and H + S goats showed significantly reduced isobutyrate and valerate concentrations compared to HGW-SARA goats (P < 0.05). Isovalerate levels were significantly decreased in HGC-Health and S + H goats relative to HGW-SARA (P < 0.01), whereas other VFAs remained comparable across all groups (Fig. 6B). Notably, LPS content was significantly higher in HGW-SARA goats compared to HGC-Health, S + H, and H + S (P < 0.01), with S + H and H + S goats showing reduced LPS levels compared to HGC-Health (P < 0.05), while NH₃-N levels showed no significant variation among the groups (Fig. 6C). Post-mortem assessments of rumen epithelial papillae indicated a significantly thicker stratum granulosum in HGW-SARA goats compared to HGC-Health, S + H, and H + S goats (P < 0.05, Fig. 6D). Clear detachment of the stratum corneum and widened gaps between layers were visually confirmed in HGW-SARA goats (Fig. 6E).

Transcriptomic profiling of rumen epithelium in HGW-SARA vs. S + H and HGW-SARA vs. H + S comparisons aimed to elucidate the mechanisms underlying the influence of the microbiome and metabolites on SARA-susceptibility and tolerance. KEGG analysis revealed four enriched immune-related pathways in HGW-SARA vs. S + H, including Natural killer mediated cytotoxicity and IL-17 signalling pathway (P < 0.001, Fig. 6F). Seven genes enriched in the IL-17 pathway, such as CCL20, IL1B, CXCL10, and CCL2 were upregulated in HGW-SARA (P < 0.05, |log2FC|> 1, Fig. 6G). At the same time, HGW-SARA vs. H + S comparison revealed 18 immune-related pathways enriched, among which Th17 Cell Differentiation, and IL-17 signalling pathway stood out (P < 0.001, Fig. 6H). Fifteen differentially expressed genes were identified in IL-17 signalling pathway, with genes TRAF5 and IL17B upregulated in H + S goats, while other pro-inflammatory genes were upregulated in HGW-SARA goats (P < 0.05, |log2FC|> 1, Fig. 6I).

The PCoA of rumen microbiome sequences among HGW-SARA, HGC-Health, S + H, and H + S goats revealed a significant difference between HGW-SARA and S + H goats (R = 0.354, P = 0.047), although no marked differences were detected in the other group comparisons (P > 0.05, Additional file 1: Fig S9A). LEfSe revealed that HGW-SARA goats harbored enriched populations of Succinivibrionaceae bacterium WG-1, Succinatimonas sp CAG 777, Ruminoccocus bromii, and Ruminoccocus sp (LDA_{HGW-SARA vs H+S} > 3, LDA_{HGW-SARA vs S+H} > 3, P < 0.05) compared to both H + S and S + H goats. Conversely, H + S and S + H goats displayed elevated levels of Bacteroidales bacterium, Prevotella sp ne3005, and unclassified-g-Prevotella (LDA_{HGW-SARA vs H+S} > 3, LDA_{HGW-SARA vs S+H} > 3, P < 0.05) relative to HGW-SARA goats (Additional file 1: Fig S9B). Functional profiling of the microbiome revealed that pathways such as Cysteine and methionine metabolism and Pyrimidine metabolism (LDA > 3, P < 0.05) to be more abundant in HGW-SARA compared to S + H goats. S + H goats, on the other hand, showed enrichment in pathways like Quorum sensing, Two-component system, Hippo signalling pathway—multiple species (LDA > 3, P < 0.05), with Two-component System also enriched in H + S goats compared to HGW-SARA (LDA = 3.21, P < 0.05, Additional file 1: Fig S9C). Comparative analysis of CAZymes via LEfSe indicated differential enrichment. For instance, CE10, CE15, and CE12 (LDA > 4, P < 0.05) were more prevalent in S + H goats, while GT92 (LDA = 5.00, P < 0.05) was dominant in HGW-SARA when comparing HGW-SARA to S + H. A similar pattern emerged with additional CAZymes when contrasting HGW-SARA with H + S goats, which found that GT92, GT2_Glyco_traf_2_4 (LDA > 4, P < 0.05) were enriched in HGW-SARA and GH64 (LDA = 4.51, P < 0.05) was enriched in H + S (Additional file 1: Fig S9D). Nine bacterial species enriched in H + S and S + H goats relative to HGW-SARA were classified across various genera and phyla, including Mycoplasma, Prevotella, unclassified-c-Clostridia, unclassified-c-Opitutae, unclassified-c-Bacteroidales, unclassified-c-Clostridiales, and unclassified-c-Verrucomicrobia, as well as the phyla Bacteroidetes, Firmicutes, Tenericutes, and Verrucomicrobia (Additional file 1: Fig S9E).

PLS-DA effectively discriminated between HGW-SARA and both S + H and H + S groups (Additional file 1: Fig S9F). KEGG enrichment analyses highlighted 12 pathways in HGW-SARA vs S + H, with 7 pertaining to metabolism, including galactose metabolism, arginine and proline metabolism, and tryptophan metabolism (P < 0.05). Similarly, 14 pathways were enriched in HGW-SARA vs H + S, 5 of which were metabolic, such as tryptophan metabolism, glycosaminoglycan biosynthesis—chondroitin sulfate/dermatan sulfate, and glycerophospholipid metabolism (P < 0.05, Additional file 1: Fig S9G). Notably, tryptophan metabolism was a common feature enriched in comparisons of HGW-SARA with both S + H and H + S goats. Analysis of specific metabolites in this pathway showed decreased 2-Formaminobenzoylacetate and increased N-Acetylserotonin and Xanthurenic Acid in H + S and S + H goats compared to HGW-SARA (Additional file 2: Table S7).

To explore the function of rumen microbiome to prevent SARA, we also constructed a SARA model of 14 dairy goats by feeding the same high RDS diets, and randomly selected 7 dairy goats to supplementarily feed B. adolescentis, which was identified as the key microbiome to help the host maintain SARA tolerance. After continuously monitoring pH in the 42nd day of the trial period, the pH of HRDS-BA dairy goats did not meet the criteria of SARA while the pH of HRDS-C dairy goats was lower than 5.8 over 3 h (P < 0.05, Fig. 7A). The milk yield of HRDS-BA goats was significantly higher than HRDS-C goats (P < 0.05, Fig. 7B). IL-6 and IL-17 were significantly decreased in HRDS-BA goats, while IL-10 increased significantly in HRDS-BA goats (P < 0.05, Fig. 7C). The concentration of total VFAs in rumen was higher in HRDS-BA goats (P < 0.05) with no other differences of ruminal VFA concentration between HRDS-C and HRDS-BA goats (Fig. 7D). However, compared with HRDS-C, the concentrations of acetate, isovalerate, valerate, and total acids in HRDS-BA goats were significantly increased (P < 0.01, Fig. 7E).

This experiment was conducted at the Animal Research and Technology Center of Northwest A&F University (Yangling, Shaanxi, China), and it was performed in accordance with the recommended guidelines from the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised 2004). The protocol was approved by the Institutional Animal Care and Use Committee at Northwest A&F University.

The rumen fluid was centrifuged at 13,000 × g for 10 min, and then used for VFAs concentration analysis. Meanwhile, the plasma samples were directly used for VFA concentration analysis. We analyzed the concentration of VFAs in rumen fluid and plasma as previously described [5]. In brief, 1 mL ruminal fluid supernatant or plasma was mixed with 200 μL metaphosphoric acid (25% w/v), after 3–4 h of standing at 5 °C, and then centrifuged for 15 min at 13,500 × g at 4 °C. Five hundred microliters supernatant was mixed with 200 µL of crotonic acid (28.45 mmol/L) and standing 0.5–2 h at 5 °C, then filtered through a 0.45-µm filter. The VFAs were separated and quantified with an Agilent 7820 A GC system equipped with a polar capillary column (AE-FFAP, 30 m × 0.25 mm × 0.33 μm) and a flame ionization detector (FID).

The quantification of free lipopolysaccharides in rumen fluid and plasma from dairy goats was carried out employing the chromogenic endpoint Tachypleus amebocyte lysate assay kit (catalog number EC80545S, Xiamen Houshiji, Ltd.) adhering meticulously to the manufacturer’s guidelines. In summary, rumen fluid samples were subjected to heat treatment at 100 °C for 30 min prior to being stored at −20 °C in preparation for lipopolysaccharide assessment. The assay was performed on a 96-well microplate with the optical density measured at 405 nm on a microplate reader (Model 3550, manufactured by Bio-Rad, Hercules, CA, USA).

The NH₃-N content of ruminal fluid was determined through the phenol-sodium hypochlorite colorimetric methodology. A standard curve was first established via a gradient series of ammonium standard solutions. In each test tube, 1.0 mL of appropriately diluted sample, standard solution, or distilled water (utilized as a blank control) was dispensed, followed by the addition of 4 mL of HCl solution and thorough mixing. A 0.2-mL aliquot of the mixture was then transferred to separate tubes, whereupon 2.5 mL each of phenol and sodium hypochlorite solutions were sequentially introduced, with subsequent blending. The tubes were subjected to a 10-min color development step in a 60 °C water bath, after which they were quenched in water, and absorbance was measured at a wavelength of 545 nm. Employing the established standard curve, the concentration of NH₄Cl in individual samples was precisely quantified. The lactic acid content in rumen fluid was assayed using the Rumen Fluid Lactic Acid Kit (A019-2-1, Nanjing Jiancheng Bioengineering Institute), adhering strictly to the kit protocol. Rumen samples, mixed with enzymatic reagents and chromogens, were incubated at 37 °C for 10 min to complete the reaction. Reaction cessation was achieved with a stopping agent, and color intensity, measured at 530 nm with a 1 cm light path, enabled precise colorimetric quantification.

Following fixation, the rumen epithelial tissues were subjected to automated dehydration, embedding, sectioning, and subsequent staining with hematoxylin for 10 to 20 min. This was succeeded by a tap water rinse lasting 1 to 3 min, differentiation in hydrochloric acid‒alcohol for 5 to 10 s, and a repeat rinse. To achieve bluening, the sections were bathed in 50 °C warm water or a gently alkaline solution until the appropriate color was obtained. After another 1 to 3 min, the sections were rinsed with tap water, incubated in 85% ethanol for 3 to 5 min and stained with matching amounts of eosin. A swift 3- to 5-s water rinse preceded a graded alcohol dehydration series, xylene clearing, and slide preparation for visualization. Microscopy images of five randomly chosen papillae were then captured at × 100 magnification using a Nikon ECLIPSE Ni-U microscope (Nikon Co. Ltd., Japan). Thereafter, with meticulous attention, the thicknesses of the stratum corneum, granular layer, spinous layer, and overall rumen epithelium within the papillae were measured and documented, contributing vital data to an exhaustive morphometric assessment.

Frozen tissue sections were triple rinsed with phosphate-buffered saline (PBS) for 5 min per cycle within a staining chamber and then subsequently submerged in a 5% film dissolving reagent for a 10-min incubation at room temperature, interspersed with another trio of 5-min PBS washes. To suppress nonspecific binding, the sections were uniformly covered with a 10% serum blocking solution and left at ambient temperature for half an hour. Thereafter, the sections were treated with primary antibodies specific for IL17-A, CD163, and KRT6A and incubated overnight at 4 °C. After a trio of 5-min washes with PBS, secondary antibodies were applied, and the samples were allowed to react at 37 °C for 30 min. Nuclear staining was accomplished using DAPI, with a subsequent 10-min incubation. This mixture was capped off with three final 5-min washes in PBS, and finally, the tablets were sealed using anti-fluorescence attenuation tablet sealant. Fluorescence images of the sections were captured using a specialized microscope camera system. Each section was first examined for clear cell visibility, and then three regions were imaged at × 400 magnification. The ImageJ software was used to analyze the images and compute the integrated density (IntDen), area, and mean gray value. The mean fluorescence intensity per sample was derived by averaging the IntDen values from the three images.

The concentrations of IL-1β, IL-6, IL-10, TNF-α, and IL-17 in plasma of HRDS-C and HRDS-BA dairy goats were detected using the respective ELISA kits (COIBO BIO, Shanghai, China) with double antibody sandwich method.

Total metagenomic DNA was extracted from rumen fluid samples using the E.Z.N.A.® Soil DNA Kit (Omega Bio-tek, Norcross, GA, USA), with repeated bead-beating following the manufacturer’s instructions. Concentration and purity of the extracted DNA were determined with TBS-380 and NanoDrop2000, respectively. The integrity of the DNA was verified through electrophoresis on a 1% agarose gel.

The extracted DNA was cut to an approximate average fragment length of 400 base pairs using the Covaris M220 (Gene Company Limited, China) for paired-end library construction. The NEXTFLEX® Rapid DNA-Seq (Bioo Scientific, Austin, TX, USA) was employed for constructing the paired-end library, wherein adapters, which were generated with sequencing primer hybridization sequences, were ligated to the blunt-ended DNA fragments. Paired-end sequencing was performed on an Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) with NovaSeq Reagent Kits according to the manufacturer’s instructions (www.illumina.com↗). After bridge PCR amplification, metagenomic sequencing was performed using an Illumina NovaSeq/HiSeq Xten (Illumina, USA) sequencing platform.

Paired-end Illumina reads underwent adaptor trimming and filtering of low-quality sequences—those shorter than 50 bp, possessing a quality score below 20, or containing ambiguous “N” bases—utilizing fastp [69] (version 0.20.0, available at https://github.com/OpenGene/fastp↗) and Trimmomatic (version 0.39, available at https://github.com/usadellab/Trimmomatic↗). Reads were aligned to the reference genome of Capra Hircus by BWA [48, 49] (version 0.7.9a, http://bio-bwa.sourceforge.net↗), and any reads mapped to the genome were discarded to avoid contamination with goat DNA. Metagenomic data assembly was accomplished with MEGAHIT [50] (version 1.1.2, found at https://github.com/voutcn/megahit↗). Contigs reaching or exceeding 300 bp in length were retained as the final assembly output, which then served as the foundation for downstream gene prediction and annotation. Metagenome sequencing generated a total of 1,785,081,216 reads, with 67,963,658 ± 1,681,281 reads (mean ± SEM) per sample. After quality control and removing host genes, a total of 1,323,991,102 reads were retained, with 67,963,658 ± 1,681,281 (mean ± SEM) per sample.

Open reading frames (ORFs) were inferred from assembled contigs utilizing MetaGene (https://github.com/CharlesJB/metagene↗). ORFs extending to or beyond 100 base pairs were extracted and translated into amino acid sequences. A nonredundant gene catalog was compiled through the application of CD-HIT [51] (version 4.6.1, https://github.com/weizhongli/cdhit/releases/tag/V4.6.1↗), which enforces a 90% sequence identity and 90% coverage threshold for clustering. Postquality control reads were mapped to the nonredundant gene catalog with 95% identity via SOAPaligner (version 2.21, https://github.com/ShujiaHuang/SOAPaligner↗), evaluating quantification of gene abundance within individual samples. The total gene abundance of each sample was normalized to the transcripts per million (TPM) value. Specifically, the gene abundance is expressed as the million-fold proportion of the gene abundance normalized by the gene length, relative to the total abundance of all genes normalized by their lengths in the sample.

For taxonomic classification, representative sequences of the nonredundant gene catalog were aligned to the NCBI NR database using Diamond [52] (version 0.8.35, http://www.diamondsearch.org/index.php↗), imposing an e-value threshold of 1e − 5. Functional annotations based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) were carried out with Diamond [52] against the KEGG database (http://www.genome.jp/kegg/↗), maintaining the same e-value cut-off. Carbohydrate-active enzymes (CAZymes) were annotated by subjecting the catalog to hmmscan (http://hmmer.janelia.org/search/hmmscan↗) analysis against the CAZy database (http://www.cazy.org/↗).

Ruminal fluid (1.5 mL) was extracted with methanol:acetonitrile (1:1) and sonicated at 40 kHz for 30 min at 5 °C, followed by protein precipitation at − 20 °C. After centrifugation, the supernatants were lyophilized and then reconstituted in 100 µL of acetonitrile:water (1:1) for UHPLC-MS/MS analysis using the UHPLC-Q Exactive HF-X system (Thermo Fisher Scientific, Massachusetts, USA) after another round of centrifugation and clarification. The LC–MS raw data were imported into the metabolomics processing software Progenesis QI (Waters Corporation, Milford, USA) for baseline filtering, peak identification, integration, retention time correction, and peak alignment. Subsequently, a data matrix of the retention time, mass‒charge ratio, and peak intensity was generated. Metabolite information was obtained by comparing the MS and MS/MS mass spectrum information with the metabolic public database HMDB.

Variance analysis was performed on the matrix file after data preprocessing. The R package ropls (version 1.6.2) was used to perform principal component analysis (PCA) and orthogonal least partial squares discriminant analysis (OPLS-DA), and 7 cycles of interactive validation were used to evaluate the stability of the model. In addition, Wilcox’ test and analysis of fold changes were performed. The selection of significantly differential metabolites was determined based on the variable importance in the projection (VIP) obtained by the OPLS-DA model and the p value of Wilcox’ test, and the metabolites with VIP > 1 and P < 0.05 were considered significantly different.

Differential metabolites among two groups were summarized and mapped into their biochemical pathways through metabolic enrichment and pathway analysis based on KEGG database. These metabolites can be classified according to the pathways they involved or the functions they performed. Enrichment analysis was usually to analyze a group of metabolites in a function node whether appears or not. The principle was that the annotation analysis of a single metabolite develops into an annotation analysis of a group of metabolites. Scipy.stats (Python package) (https://docs.scipy.org/doc/scipy/↗) was exploited to identify statistically significantly enriched pathways via Fisher’s exact test.

Total RNA was isolated from rumen epithelial tissue specimens employing TRIzol® Reagent, adhering meticulously to the manufacturer’s guidelines (Invitrogen). Genomic DNA was eliminated through treatment with DNase I (TaKara). The integrity of RNA and the absence of degradation or contaminant bands were verified via 1% agarose gel electrophoresis. Thereafter, the quality of the RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies) and quantified with a NanoDrop ND-2000 spectrophotometer (NanoDrop Technologies). Solely RNA samples meeting the premium quality criteria (OD260/280 = 1.8 ~ 2.2, OD260/230 ≥ 2.0, RIN ≥ 8.0, 28S:18S ≥ 1.0, > 1 μg) were advanced for the construction of sequencing libraries.

RNA purification, reverse transcription, library construction, and sequencing were executed in accordance with the manufacturer’s protocols (Illumina, San Diego, CA). Transcriptome library preparation was accomplished utilizing the TruSeq™ RNA Sample Preparation Kit from Illumina. Briefly, messenger RNA isolation was achieved through polyA selection employing oligo(dT) beads preceding fragmentation by a dedicated buffer. Thereafter, double-stranded cDNA was generated using a SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen, CA) in conjunction with random hexamer primers (Illumina). The resulting cDNA underwent end repair, phosphorylation, and adenylation according to Illumina’s recommended library preparation steps. The size of the 300 bp cDNA fragments was selected on a 2% low-range Ultra agarose gel, followed by 15 cycles of PCR amplification with Phusion DNA Polymerase (NEB). Libraries were quantitated by TBS380, and the paired-end RNA-seq library was subsequently sequenced with an Illumina NovaSeq 6000 sequencer (2 × 150 bp read length).

The raw paired-end reads were trimmed and subjected to quality control using fastp [69]. Next, the cleaned reads were mapped to the reference genome utilizing HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtmL↗) [70] for orientation alignment. The mapped reads from each sample were assembled with reference guidance by means of StringTie (https://ccb.jhu.edu/software/stringtie/↗) [71].

To discern differential expression genes (DEGs) between distinct groups, gene expression levels were quantified employing the transcripts per million reads (TPM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/↗) [72] was employed for gene expression abundance estimation. Differential expression analysis was primarily conducted with DESeq2 [73], where DEGs with a |log2 (fold change)|≥ 1 and a significance level of P value ≤ 0.05 were deemed significantly differentially expressed. Moreover, enrichment analyses for Gene Ontology (GO) terms and KEGG pathways were conducted to pinpoint DEGs showing significant enrichment (P value ≤ 0.05). This functional enrichment utilized Goatools (https://github.com/tanghaibao/Goatools↗) [74] and KOBAS (http://kobas.cbi.pku.edu.cn/home.do↗) [75] for GO and KEGG analysis.

After rinsing with PBS, the rumen epithelial tissue samples stored in liquid nitrogen were subjected to nucleus extraction to prepare single-cell nuclear suspensions. The 10x Genomics ChromiumTM system was employed to generate Gel Bead-in-Emulsion (GEMs), with subsequent GEM collection for reverse transcription and barcoding within a PCR machine. Following lysis of the GEMs, first-strand cDNA was enriched using magnetic beads and purified, followed by cDNA amplification and quality control. Qualifying cDNA underwent next-generation sequencing library preparation, involving fragmentation, adapter ligation, indexing PCR, and final library quantification and quality checks. The paired-end 150 bp (PE150) qualified libraries were subsequently sequenced on the Illumina HiSeq platform.

The quality of the raw read data was analyzed by Fastp [69]. Cell Ranger was used for quality control of the raw data, and the comparison software STAR [76] was used to compare the sequenced reads to the reference genome of Capra Hircus to obtain quality control results, such as the number of high-quality cells, number of genes, and genome comparison information of the sample data. Seruat was used for further quality control filtering and removal of outliers on the basis of initial quality control and standardization in “LogNormalize” [77, 78]. After dimensionality reduction using UMAP, the “SNN” method in Seruat was used for clustering. SingleR was used for cell type identification, and epithelial cell subclusters were identified with reference to the bovine single-cell database of Zhejiang University (http://cattlecelllandscape.zju.edu.cn↗). All the visualization of UMAP diagrams were based on R packages of Seruat, ggplot2, ggunchull, ggrepel, and tidyverse.

An appropriate quantity of rumen epithelial tissue was harvested and thoroughly homogenized before lysis. Total RNA extraction was performed utilizing TRIzol reagent (Carlsbad, USA). Thereafter, the extracted RNA was reverse transcribed into cDNA following the standard protocol provided by the PrimeScript® RT Kit (Takara, China). Quantitative real-time polymerase chain reaction (qRT‒PCR) assays were executed with SYBR® Premix Pro TaqTM II (Takara, China) on a real-time fluorescence quantitative thermal cycler (Light Cycler 9,603,030,973). Each 10 μL qRT‒PCR mixture comprised 4 μL of cDNA template, 0.5 μL each of forward and reverse primers, and 5 μL of SYBR Green mix. The PCR cycling conditions entailed an initial denaturation step at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for annealing and extension. The specific primer pairs utilized for the target genes are enumerated in Additional file 2: Table S10. Gene expression levels were analyzed employing the widely accepted approach for relative quantification in qRT‒PCR experiments known as the 2^−ΔΔCt method [79].

The structure equation model was constructed in R studio with the R package lavaan [80] to explore the relationships among the rumen microbiome, rumen metabolome, gene expression of the rumen epithelium, and number of epithelial cells. A variance‒covariance matrix was constructed, and the model parameters were estimated utilizing the maximum likelihood method. Prior to commencing the process, the abundance data underwent a transformation to adhere to a normal distribution, specifically adding one to the original value and taking the logarithm. According to the model’s fitting effectiveness, which is referred to as the comparative fit index (CFI), we iterated the model construction several times and compared it with ANOVA. After building successful models, view standardized path coefficients and visualize them.

The power analysis was performed using the R packages "samr" and "limma" on integrated phenotypes alteration, RNA-seq dataset metagenomic datasets and metabolomic datasets to determine the optimal sample size [81]. The results showed that 5 biological replicates per group were sufficient to achieve a statistical power of 0.831, meeting the standard threshold (power > 0.8). The distance algorithm was used to analyze the species and functional composition among samples by PCoA, and the differences between groups were tested by ANOSIM. LEfSe was used to analyze differential microorganisms and their functions, and ggplot2 was used for visualization. DESeq2 software was used to screen DEGs, and KEGG enrichment analysis of DEGs was performed by KOBAS. Differentially abundant metabolites were determined based on the variable weight value (VIP) obtained by the OPLS-DA model and the P value of Wilcox’ test, and statistical analysis of differentially abundant metabolites was performed based on Wilcox’ test.

Apart from above, all the statistical analyses were conducted in IBM SPSS Statistics 27 using the mean ± SEM and one-way ANOVA followed by the LSD and Duncan tests. *P < 0.05, **P < 0.01, ***P < 0.001 indicate significance and were visualized with GraphPad Prism 9.5. The correlation analysis was performed in R studio with R packages such as BiocManager, ComplexHeatmap, and circlized with Pearson or Spearman correlation coefficients. ROC curve analysis was also conducted in GraphPad Prism 9.5 with the Wilson/Brown method and 95% confidence intervals.

Andrew Cosgrove was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team. The peer-review history is available in the online version of this article.

This experiment was conducted at the Animal Research and Technology Centre of Northwest A&F University (Yangling, Shaanxi, China), and it was performed in accordance with the recommended guidelines from the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised 2004) with approval number: NWAFU-DK2022209.

Not applicable.

The authors declare no competing interests.

All the data generated or analysed for this study are included in this paper. The metagenome data, RNA-seq data, metabolomics data and single-nucleus RNA-seq data were deposited into the China National Center for Bioinformation (CNCB; https://www.cncb.ac.cn/?lang=en↗) under accession number PRJCA026627 [82].