Transcriptional adaptation of rumen papillae to high-grain diet reveals distinct temporal phases and SARA susceptibility signatures

TL;DR

This study explores how the rumen adapts to high-grain diets and identifies genes linked to susceptibility to subacute ruminal acidosis in cows.

Contribution

The study reveals distinct temporal phases of transcriptional adaptation and identifies potential biomarkers for SARA susceptibility in cattle.

Findings

Rumen papillae show biphasic transcriptional adaptation to high-grain diets, with early and late response phases.

Genes like CCDC196 and MYO7B are under-expressed in SARA-susceptible cows, suggesting roles in acidosis predisposition.

SARA-resistant cows exhibit stronger transcriptome-metabolome correlations, indicating more coordinated adaptation.

Abstract

Transitioning to a high-grain (HG) diet significantly alters rumen fermentation by increasing the production of short-chain fatty acids (SCFAs) and lowering rumen pH, which may contribute to subacute ruminal acidosis (SARA) and damage to the ruminal epithelium. Rapid adaptation of rumen epithelium to these metabolic shifts is essential to maintain homeostasis, but the transcriptional mechanisms underlying this adaptation remain poorly understood.

We analyzed the temporal progression of gene expression and metabolomic profile in rumen papillae collected during low-grain feeding (LG) and one week after transitioning to a HG diet (HG1), or four weeks after (HG4) in cows classified as susceptible or resistant to SARA. RNA sequencing identified 955 differentially expressed genes (DEGs) across time points, revealing a biphasic adaptation pattern. Early responses (HG1) showed moderate transcriptional changes, while HG4 was characterized by substantial transcriptional remodeling. Pathway analysis indicated three major functional categories affected during adaptation: cellular stress response, metabolic adaptation, and protein processing. Notably, sterol biosynthesis genes showed transient upregulation at HG1 followed by downregulation at HG4, coinciding with morphological changes in rumen wall thickness and n-butyrate concentration in rumen fluid. Correlation analyses comparing gene expression patterns and metabolite level changes triggered by the dietary transition revealed potential links between metabolic and transcriptional adaptation. Of particular interest, valerate levels at HG1 correlated with genes involved in tissue remodeling at HG4, implying that valerate may contribute to delayed epithelial responses. Next, transcriptional differences between SARA-susceptible and SARA-resistant animals included genes related to inflammation, cell structure, and metabolism that persisted across all time points, suggesting underlying intrinsic differences in SARA susceptibility that are present before and persist during dietary challenge. Key genes consistently differentially under-expressed in SARA-susceptible animals, CCDC196 and MYO7B, represent potential biomarkers for SARA predisposition. Finally, the SARA-resistant group showed a greater number of transcriptome-metabolome correlations, suggesting more coordinated epithelial responses to diet change compared to the SARA-susceptible group.

Our findings provide insights into the molecular mechanisms underlying rumen adaptation to HG diets and individual variation in SARA susceptibility, providing a basis for developing strategies to optimize dietary transitions in ruminant production systems.

The online version contains supplementary material available at 10.1186/s40104-026-01352-8.

Figures

Temporal gene expression changes in rumen papillae during high-grain diet adaptation. A Experimental timeline. B Cows were transitioned from a low-grain (LG, 40% concentrate) diet to a high-grain (HG, 65% concentrate) diet, with tissue collected at LG, HG1 (1 week), and HG4 (4 week) time points. Heatmap of the top 50 differentially expressed genes (DEGs) based on adjusted *P*-value across LG, HG1, and HG4. Hierarchical clustering was performed using ward.D2 and euclidean distance. C The overlap of DEGs between pairwise comparisons (LG vs. HG1, LG vs. HG4, HG1 vs. HG4). DEGs were defined by an

Functional pathways affected during rumen adaptation to a high-grain diet. A Pathway enrichment analysis for DEGs across time points was performed using the KEGG database. The top 20 pathways ranked by adjusted *P*-value are shown. B The resulting pathways were clustered into three categories: metabolic adaptation, cellular stress response, and protein processing. DEGs assigned to KEGG pathways were extracted, and their fold expression changes (log_2_FC) were compared among the three time points. C The changes were visualized with the category corresponding to each DEG. Expression patterns of

Metabolic changes associated with rumen adaptation. A–M Ruminal short-chain fatty acid (SCFA; A–D) as well as other indicated metabolite concentrations (E–M) across the selected time points. N Biplot showing the first canonical variate for CCA analysis of merged metabolic and transcriptomic data. O–P Gene-metabolite correlation network (|*r*| > 0.7, FDR < 0.05) for time point HG1 (O) and HG4 (P). R Red lines indicate negative and green lines positive correlation. The results of gene ontology “biological process” analysis for genes correlating with valerate (|*r*| > 0.7 and FDR < 0.05)

Structural changes associated with rumen adaptation. A Representative histological sections of rumen papillae at each time point (H&E stain). B–D Quantification of epithelial parameters: stratum corneum thickness, epithelial thickness, and rumen wall thickness (RWT). E Normalized count for *IGFBP3* expression across time points. Data are presented as mean ± SD. ^*^*P* < 0.05, *n* = 10 or 11. Depending on the data distribution, ANOVA with Tukey’s or Kruskal–Wallis with Dunn’s post-hoc tests were applied to verify statistical significance

Transcriptomic signatures of SARA susceptibility in rumen papillae. A Heatmap of 96 DEGs between SARA-susceptible (S) and SARA-resistant (R) cows across all time points (*P* < 0.01, log_2_FC > ± 1.5). B PCA trajectory plot of SARA S and R cows across three time points. Individual trajectories (thin lines) and group means (bold lines) are shown

Functional categorization of genes differing between SARA phenotypes. Bar graph summarizing pathway enrichment analysis of DEGs between the SARA S and SARA R groups, classified into inflammation, metabolism (nutrient transport, lipid, and protein metabolism), and cell cycle/cell structure. The shape of the figure indicates the time point, and its size reflects the *P*-value for each time point

Inflammatory and metabolic differences between SARA phenotypes. A STRING protein–protein interaction network for DEGs related to inflammation between SARA groups. B SCFA concentrations in rumen fluid from SARA S and R animals at the indicated time point. C–F The indicated metabolite concentrations at LG and HG1 time points. Data are expressed as mean ± SD; ^*^*P* < 0.05

Correlation of metabolome and transcriptome in SARA groups. A Volcano plot indicating Spearman Rho and *P*-value for correlation of genes and metabolites for R and S groups. B Correlations specific to R are on the left and to S on the right side of the graph. Chord diagram visualizing gene-metabolite correlation with |Δρ|> 0.5 and *P* < 0.05 threshold. C Red lines indicate connections for group S and green for group R. Pathway enrichment analysis using WIKI pathways comparing DEGs correlating with metabolites (|Δρ| > 0.5 and *P* < 0.05) in R and S groups

Full text

Introduction

Ruminant production globally is increasingly adopting high-grain (HG) diets to meet the growing demand for animal products. However, the transition from forage-based to grain-rich diets poses significant challenges to rumen health and function due to increased fermentation, alterations in short-chain fatty acid (SCFA) profiles, and a reduction in ruminal pH [1]. The diet-induced changes can lead to subacute ruminal acidosis (SARA), a prevalent digestive disorder characterized by periods when ruminal pH falls below 5.6–5.8 for several hours each day [2, 3]. SARA induces structural damage to the ruminal epithelium, promoting papillary erosion and inflammatory infiltration, which compromises epithelial barrier integrity and facilitates bacterial translocation, particularly of lipopolysaccharide endotoxins, into the systemic circulation [4, 5]. Subsequently, it initiates a cascade of inflammatory responses and contributes to liver damage, including abscesses, through the activation of hepatic immune responses [4, 6]. These systemic effects mediated by ruminal dysbiosis, endotoxemia, and inflammatory responses extend beyond the gastrointestinal tract and result in reduced feed intake, milk yield depression, laminitis, diarrhea, impaired immunity, and decreased reproductive performance [4, 6, 7].

The rumen epithelium acts as a selectively permeable barrier, regulating nutrient absorption while also protecting the underlying tissues from potentially harmful substances in the lumen [8, 9]. To maintain these functions and preserve homeostasis, the rumen epithelium must quickly adapt to the diet changes. The epithelium coordinates responses at morphological, functional, and molecular levels due to its remarkable adjustment capacity. HG diets generally lead to the growth of papillae, which increases the absorptive surface area to accommodate the absorption of the increasingly available SCFAs [10, 11]. At the cellular level, the adaptation of the rumen epithelium involves the modulation of tight junction proteins, an increase in SCFA transport capacity, and enhanced metabolic activity [12].

Although significant progress has been made in understanding rumen adaptation, the temporal dynamics of the transcriptional responses involved in this process remain incompletely characterized. Previous studies have identified important changes in gene expression related to metabolic functions, cell proliferation, and immune responses during adaptation to HG diets [13–15]. Some have suggested that this adaptation occurs in distinct phases, with immediate responses aimed at maintaining pH homeostasis, followed by longer-term structural and functional modifications [5]. However, the regulatory networks controlling these responses, especially extending beyond the third week of the diet, have not been fully elucidated.

The main objective of this study was to characterize the transcriptional dynamics of rumen papillae as they adapt to an increased dietary grain content. To enhance mechanistic understanding, the correlations between transcriptome and rumen fluid metabolite data were assessed. We analyzed samples collected at three different points: before the dietary transition (low grain, LG), one week after the transition (high grain 1, HG1), and four weeks after the transition (high grain 4, HG4). In addition to demonstrating time-dependent adaptation to an HG diet, we aimed to identify transcriptional signatures that distinguish SARA-susceptible animals. This investigation aimed to provide novel insights into the molecular basis of rumen adaptation and to inform the development of strategies to optimize dietary transitions in production settings.

Materials and methods

Experiment design and sample collection

The experimental procedure was approved by the Ethics and Animal Welfare Committee of the University of Veterinary Medicine Vienna and the Austrian Ministry for Education, Science, and Research (protocol number: BMBWF-2022-0.276.659). The trial was conducted using eleven Holstein cows in their early-mid second lactation (646 ± 59 kg of average body weight and 87 ± 57 d in milk at the start of the experiment, body condition scoring 3.0 ± 0.2). The animals were rumen-cannulated and housed in a free‐stall barn at the research farm (Pottenstein, Austria) of the University of Veterinary Medicine, Vienna, between July and November 2022. The animals were fed a control diet with 40% grain feeding (low grain; LG) (Fig. 1A). The cows were transitioned to a HG (65%-grain) diet for the following four weeks. The diet composition is presented in Table S1.Fig. 1. Temporal gene expression changes in rumen papillae during high-grain diet adaptation. A Experimental timeline. B Cows were transitioned from a low-grain (LG, 40% concentrate) diet to a high-grain (HG, 65% concentrate) diet, with tissue collected at LG, HG1 (1 week), and HG4 (4 week) time points. Heatmap of the top 50 differentially expressed genes (DEGs) based on adjusted P-value across LG, HG1, and HG4. Hierarchical clustering was performed using ward.D2 and euclidean distance. C The overlap of DEGs between pairwise comparisons (LG vs. HG1, LG vs. HG4, HG1 vs. HG4). DEGs were defined by an adjusted P-value < 0.05 (DESeq2)

The animals were classified into SARA-susceptible and SARA-resistant groups based on rumen pH and individual SARA severity during their first lactation as described earlier by Hartinger et al*.* [16]. The SARA-susceptible group, consisting of five animals, experienced a ruminal pH < 5.8 for 45% of the time during their first lactation, and the SARA-resistant group, consisting of six animals, experienced a ruminal pH < 5.8 for 8% of the time during their first lactation [16]. Data for mean, maximum, and minimum pH, as well as time with mean pH below 5.8 per day, are presented in Supplementary Fig. S1A.

To determine the sample size, a power analysis was conducted using data from the previous experiment for nine significantly altered proteins (1.4–4.4 fold changes; P = 0.0002–0.0247). The measured statistical power for these effects at n = 6 per group ranged from 70% to 99.9%. Based on these data, n = 6 animals per group was selected as it provides adequate power to detect biologically relevant changes (≥ 1.5-fold) while following 3Rs principles for animal research.

Rumen papillae samples were collected at baseline during LG feeding, one week (HG1), and four weeks (HG4) after the transition to HG was completed. Free rumen liquid samples were gathered at time point LG and each week of HG feeding (HG1, HG2, HG3, and HG4). Both types of samples were collected 4 h after the morning feeding using the methods described previously [17]. Prior to the collection of rumen papillae biopsies, the rumen was partially emptied. The tissue samples for gene expression analyses were snap-frozen and then stored at −80 °C.

Histological samples were fixed in a formalin solution (4% v/v, ROTI^®^Histofix, Carl Roth, Germany) for 24 h and stored in 70% ethanol until further processing.

Immediately after sampling, the rumen was refilled with its contents. Each part of the equipment was washed and thoroughly disinfected after every use.

Rumen fluid for metabolomics was collected weekly from the ventral sac of the rumen using a sterile 20-mL syringe 4 h after morning feeding. All samples of rumen content were snap-frozen and stored at −80 °C.

Blood was collected weekly from the jugular vein using 9 mL vacutainer tubes (Vacuette, Greiner Bio-one, Kremsmünster, Austria) before offering fresh feed in the morning. Blood samples were centrifuged at 2,000 × g for 10 min at 4 °C after clotting at room temperature for 1.5 h, and the collected serum was stored at −80 °C until further analysis.

RNA extraction and sequencing

RNA extraction was conducted as previously described [18]. Briefly, the RNeasy Mini QIAcube Kit (Qiagen, Germany) was used with minor protocol modifications. Approximately 25 mg of papillae was combined with 350 µL of RLT buffer in 2 mL tubes containing 0.6 g of ceramic beads. Using a Fastprep-24 instrument (MP Biomedicals, USA), the samples were homogenized and then centrifuged at 10,000 × g for 1 min. The lysate was transferred to a 2-mL tube and centrifuged again at 14,680 × g for 3 min. Subsequent steps were executed as specified by the manufacturers, though centrifugation durations were extended from 15 to 30 s. After adding 500 µL of buffer RPE, the column was centrifuged at maximum speed for 1 min to dry the membrane. The filter was then placed into a new 1.5-mL tube and allowed to dry for 1 min. In the final step, 30 µL of RNase-free water was added to the filter and incubated for 1 min. RNA was eluted by centrifugation at 10,000 × g for 1 min and stored at −20 °C. Genomic DNA was eliminated using DNAse I (Ambion^®^ TURBO DNA-free), after which RNA integrity was evaluated using the Qubit RNA IQ Assay Kit, and the extracted RNA was quantified with the Qubit RNA HS Assay Kit (Invitrogen, Thermo Fisher Scientific, Austria) on the Qubit Fluorometer 4.0 (Life Technologies Corporation, USA).

The following steps were outsourced and performed by CeGaT GmbH (Tübingen, Germany). Total RNA (100 ng/sample) was used for library preparation using the KAPA RNA HyperPrep Kit with RiboErase (HMR) in combination with KAPA Globin Depletion Hybridization Oligos (Roche, Austria). This protocol depletes ribosomal and globin RNA, enabling enrichment of coding and non-globin transcripts. Libraries were sequenced on an Illumina NovaSeq 6000 platform using a paired-end strategy with read lengths of 2 × 101 bp. The additional base in both reads was sequenced to improve the quality score accuracy for the terminal nucleotide. Sequencing quality was high, with Q30 values of ≥ 87.83%. Demultiplexing of the sequencing reads was performed using Illumina bcl2fastq (version 2.20). Adapters were trimmed with Skewer (version 0.2.2) [19]. Quality trimming of the reads has not been performed. Trimmed raw reads were aligned to Bos taurus ARS-UCD1.2 using STAR (version 2.7.3) [19].

Sequencing data analysis

Principal components analysis (PCA) was performed to identify the major factors responsible for distribution variation. The analysis was conducted using the ggplot2 and ggrepel packages in R (version 4.0.4; R Core Team 2015). Additionally, Cook's Distance was calculated to assess outliers. Upon data visualization, one sample from the LG time point was identified as an outlier (Fig. S1B), and the Cook's Distance test confirmed this. The rest of the analysis was performed excluding this sample.

Differential expression analysis between groups was performed with DESeq2 (version 1.24.0) [20]. DESeq2 uses a negative binomial generalized linear model to test for differential expression based on gene counts. An overall Likelihood Ratio Test (LRT) was performed to detect genes with significant expression changes across all groups. Log_2_ fold changes (log_2_FC), P-values, and adjusted P-values (Benjamini–Hochberg false discovery rate, FDR) were calculated for each comparison. For further analyses, only differentially expressed genes (DEGs) with an adjusted P-value < 0.05 were considered statistically significant. Next, splitting the three diet groups into six diet-SARA groups reduced the number of biological replicates, reducing statistical power. For that reason, cutoff criteria in the analysis of the differences in RNA expression were reduced to P-value < 0.01 and |log_2_FC| > 1.5.

A post hoc power analysis was conducted to assess the sensitivity of the RNA-seq analysis. Based on the sample sizes and the effect size distribution across over 15,000 genes, the study had 23.1% power for LG vs. HG4 comparisons, 22.8% for HG1 vs. HG4, and 13.8% for LG vs. HG1 at α = 0.05 (Fig. S1C). This analysis revealed that while the study design was sufficient to detect medium-to-large differential expression effects (mean |Cohen's d| = 0.41–0.58), the statistical power was limited for detecting smaller effects.

To test for differences in community composition among groups, we performed permutational multivariate analysis of variance (PERMANOVA) using the adonis2 function and analysis of similarities (ANOSIM) in the vegan package in R. Both tests were conducted with 999 permutations based on Bray–Curtis dissimilarity matrices. Statistical significance was assessed at α = 0.05.

Hierarchical clustering analysis was conducted to visualize the relative expression patterns of genes across different experimental groups. The analysis was performed using the pheatmap and RColorBrewer packages, applying the ward.D2 clustering method. To evaluate the clustering quality, the Cluster and Stats packages were used to calculate the cophenetic correlation and average silhouette width.

Functional enrichment analysis of DEGs was conducted using ShinyGO (version 0.80; http://bioinformatics.sdstate.edu/go80/) against the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database and the Gene Ontology (GO) Biological Process ontology. The analysis was performed using the default parameters with FDR correction. Pathways with FDR < 0.05 were considered significantly enriched. The protein interaction network was visualized using the STRING database (https://string-db.org/).

Histology

Biopsies of the ruminal villi were collected via rumen fistulas as described in Section “Experiment design and sample collection”. The samples were fixed in formalin for 24 h and subsequently stored in 70% ethanol at 4 °C. The papillae were then dehydrated, embedded in paraffin, and sectioned into 5 µm-thick slices. For hematoxylin and eosin staining (hematoxylin–eosin stain, Carl Roth, Germany), the slides were deparaffinized and stained. Images were captured using a Zeiss Axioskop and an Axiocam 504 (Carl Zeiss, Jena, Germany) with ZEN 2.3 lite software (Zeiss ZEN^®^, Carl Zeiss AG).

Metabolomics analyses

SCFA concentrations in rumen liquid samples from LG, HG1, HG2, HG3, and HG4 time points were determined using gas chromatographic analysis following established protocols [21]. Samples were centrifuged (20,000 × g, 25 min, 4 °C) to remove residual digesta particles. Afterwards, 0.6 mL of the clarified supernatant was collected into clean tubes. Sample preparation involved acidification with 0.2 mL of 1.8 mol/L HCl, followed by the addition of 0.2 mL of internal standard solution (4-methylvaleric acid, Sigma-Aldrich, St. Louis, MI, USA). Following a second centrifugation step (20,000 × g, 25 min, 4 °C), the resulting supernatant was collected in GC vials for chromatographic analysis. SCFA quantification was conducted using a Shimadzu GC Plus system equipped with flame-ionization detection and a 30 m × 0.53 mm ID × 0.50 μm column (Trace TR Wax, Thermo Fisher Scientific).

Carboxylic acids, sugar phosphates, and nucleotides were determined in LG and HG1 rumen liquid samples by anion-exchange chromatography-high resolution mass spectrometry (IC-HRMS) on a Dionex Integrion HPIC system (Thermo Fisher Scientific) coupled to a Thermo Fisher Scientific Q Exactive Orbitrap mass spectrometer as described in detail by Ricci et al*.* [22]. In short, sample preparation consisted of shaking 20 µL of rumen fluid with 980 µL of acetonitrile/water (80:20, v/v) at 4 °C for 10 min, centrifugation and the supernatants were diluted 1:10 with acetonitrile/water (80:20, v/v). Analytes were separated at 30 °C on a Dionex IonPac AS11-HC column (250 mm × 2 mm, 4 μm) with a matching AG11-HC guard column (50 mm × 2 mm, 4 μm) using a KOH gradient. Potassium ions were then exchanged for protons using a Dionex AERS 500 suppressor (2 mm) operated at 76 mA. Mass spectrometric detection was carried out in full scan mode (m/z 50–750) at a resolution of 70,000 after electrospray ionization (ESI) in negative mode.

The other analytes were determined in rumen liquid samples from LG and HG1 time points by high-performance liquid chromatography (HPLC, one reversed-phase (RP) method and one hydrophilic interaction chromatography (HILIC) method) coupled to tandem mass spectrometry (MS/MS) as described in detail by Xu et al. [23]. Rumen fluid was worked-up by vortexing 200 µL of rumen fluid with 800 µL of cold acetonitrile for 10 s and centrifugation at 14,350 × g at 4 °C for 10 min. Supernatants were further diluted 1:20 (v:v) with acetonitrile/water (80:20, v:v). The undiluted supernatant and the additional dilution were analyzed using both HPLC–MS/MS methods. HPLC–MS/MS analyses were carried out on an Agilent 1290 Infinity II series UHPLC system (Agilent, Waldbronn, Germany) coupled to a 6500+ QTrap mass spectrometer equipped with an IonDrive Turbo V^®^ source (Sciex, Foster City, CA, USA). RP chromatography was performed on a Kinetex C18 column (150 mm × 2.1 mm, 2.6 µm, Phenomenex, Aschaffenburg, Germany) in gradient elution mode. HILIC separations were carried out on a Kinetex HILIC column (150 mm × 2.1 mm, 2.6 µm, Phenomenex), also with gradient elution. Mass spectrometric detection was done in scheduled selected reaction monitoring mode with fast polarity switching after electrospray ionization.

Analytes were quantified on the basis of calibration functions measured at the beginning and end of each measurement sequence, along with regular injections of a pooled quality control (QC) sample and a single calibration standard throughout the run. Drift correction was performed using loess regression on the QC sample or standard and peak areas were normalized accordingly. Subsequently, quantification was performed using linear or quadratic regression models based on manual revision. All steps were carried out with the QuantyFey application [24].

Metabolomics and transcriptomics data correlation

To explore relationships between transcriptomic and metabolomic data, canonical correlation analysis (CCA) was performed. RNA-seq and metabolite matrices were scaled and dimensionally reduced using PCA, retaining components that explained ≥ 80% of the variance. CCA was then applied to the reduced datasets to extract canonical correlations, which were evaluated for statistical significance using Wilks’ Lambda. Correlation results were visualized via scatterplot of canonical variates.

Spearman correlation analysis was performed between all differentially expressed genes (995 genes) and quantified metabolites (66 metabolites) using matched samples from the same animals. To address multiple testing while maintaining sensitivity for biologically relevant associations, a two-stage FDR correction approach was employed. First, we identified all correlations with absolute Spearman's rho > 0.5, representing moderate to strong effect sizes. Second, Benjamini–Hochberg FDR correction was applied only to this pre-filtered set of biologically plausible correlations. This approach reduces the multiple testing burden while focusing on correlations most likely to represent meaningful biological relationships. Correlation significance was defined as FDR < 0.05 in the two-stage analysis. Gene-metabolite pairs with |r| > 0.7 were considered strong correlations and presented in the network figures.

For SARA groups, differential correlation analysis was conducted between transcriptomic and metabolomic profiles in SARA-susceptible (S) vs. SARA-resistant (R) groups, identifying significant metabolite-gene pairs based on changes in Spearman’s rank correlation (|Δρ| > 0.5, P < 0.05). A volcano plot was generated to highlight the top differential correlations and candidate biomarker relationships. To visualize high-confidence correlations between metabolites and genes, a chord diagram was generated by filtering a list of metabolite-gene pairs with |Δρ| > 0.5 from a differential correlation analysis.

Data processing and visualization were performed using the readxl, dplyr, tidyr, ggpubr, reshape2, ggrepel, ggplot2, circlize, and tibble packages.

Statistical analyses

For histology data and metabolite levels, statistical analyses were conducted using R, primarily with the tidyverse, rstatix, and ggpubr packages. For comparing two groups, Student’s t-test was applied. For data sets containing three groups, the normality of data distribution within each group was assessed using the Shapiro–Wilk test. For datasets with normally distributed values across groups, one-way analysis of variance (ANOVA) was conducted, followed by Tukey's Honest Significant Difference (HSD) test for post-hoc pairwise comparisons. For non-normally distributed data, the Kruskal–Wallis test was used, followed by Dunn’s test for multiple pairwise comparisons with Bonferroni adjustment. Outliers were identified using the interquartile range (IQR) method. A statistically significant difference was considered to be present at P < 0.05.

Results

Rumen papillae exhibit a biphasic transcriptional response to high-grain (HG) diet

PCA analysis revealed modest sample separation along the principal components, with PC1 explaining 19.4% and PC2 explaining 15% of the total variance (Fig. S1D). The distribution pattern suggests that PC1 primarily pictures the temporal progression of adaptation, with HG4 samples generally separated from LG samples along this axis. PC2 reflected individual animal variability in response to dietary change, hinting at adaptation heterogeneity. Multivariate community composition differed significantly among the three groups (PERMANOVA: F₂ = 2.362, R^2^ = 0.14, P = 0.022). Analysis of similarities confirmed significant differences in community structure between groups (ANOSIM: R = 0.153, P = 0.012).

Analysis using DESeq2 resulted in 995 statistically significantly affected (adjusted P-value < 0.05) DEGs across all time points. Hierarchical clustering analysis was conducted to visualize the relative expression patterns of genes across different experimental groups. Clustering the top 50 most significantly differentially expressed mRNAs revealed distinct expression patterns across time points, with an average silhouette width of 0.34 and a cophenetic correlation coefficient of 0.94 (Fig. 1B). The dendrogram structure for samples demonstrated a clear separation between HG1 and HG4 groups, indicating substantial transcriptional differences between early and late adaptive responses. LG samples showed closer clustering with HG1 time point (Fig. 1B). The side dendrogram was divided into two major gene clusters with distinct temporal expression patterns: an upper cluster showing primarily upregulation in HG4 and a lower cluster picturing downregulation in HG4 compared to LG and HG1.

Pairwise comparison of DEGs between time points revealed substantial variations across the three time points (Fig. 1C, Table S2). The comparison between HG1 and HG4 showed the highest number of DEGs (697 DEGs), representing the majority of differentially expressed genes, followed by LG vs. HG4 (422 DEGs), and relatively few transcriptional differences in LG vs. HG1 (22 DEGs). These temporal patterns of DEG accumulation suggest distinct biological phases of rumen epithelial adaptation. The extensive reprogramming between HG1 and HG4 indicates a delayed adaptive response upon sustained metabolic pressure. This contrasts with the limited transcriptional changes between LG and HG1.

Interestingly, only one gene (ADPRM) was found to be significantly regulated in both the LG vs. HG1 and LG vs. HG4 comparisons (Fig. 1C, Table S2), suggesting very distinct transcriptional responses between short- and long-term adaptation to diet change. In both cases, ADPRM was upregulated compared to LG. Among the 10 genes common to LG vs. HG1 and HG1 vs. HG4 comparisons (CDHR4,* C7*,* PLPP1*,* ACSBG1*,* PARP3*,* BPIFA1*,* NGEF*,* INKA1*,* BMP6*, and MST1R), all were regulated in opposite directions (Table S2). Opposite regulation of genes shared across multiple comparisons suggests that pathways are dynamically recalibrating as the system shifts from an acute response to long-term adaptation. Of the 199 genes simultaneously differentially expressed in LG vs. HG4 and HG1 vs. HG4 comparisons, 122 showed consistent directional changes, while 77 were regulated in opposite directions. Overall, the analysis showed the strongest contrast in DEG expression between HG1 and HG4 time points.

Dietary adaptation primarily affects metabolic pathways, stress response, and protein processing

The 995 differentially expressed transcripts were matched against the KEGG database to assess which processes are involved in adaptation to the diet change. Figure 2A and corresponding Table S3 present some of the most statistically significantly affected pathways.Fig. 2. Functional pathways affected during rumen adaptation to a high-grain diet. A Pathway enrichment analysis for DEGs across time points was performed using the KEGG database. The top 20 pathways ranked by adjusted P-value are shown. B The resulting pathways were clustered into three categories: metabolic adaptation, cellular stress response, and protein processing. DEGs assigned to KEGG pathways were extracted, and their fold expression changes (log_2_FC) were compared among the three time points. C The changes were visualized with the category corresponding to each DEG. Expression patterns of DEGs involved in lipid metabolism and steroid/terpenoid backbone biosynthesis

The mRNAs from the "steroid biosynthesis" term emerged as a particularly prominent pathway, characterized by a high odds ratio (Fig. 2A). Terpenoid backbone biosynthesis follows closely, suggesting a further contribution to sterol metabolism.

Nerve cell-related pathways, including “Pathways of neurodegeneration—multiple diseases”, “Herpes simplex virus 1 infection”, and “Axon guidance” pathways, showed significant enrichment (Fig. 2A). Importantly, DEGs assigned to these pathways included genes that may be expressed by neurons, but none of them are neuron-specific, including SEMA3B/4C/6a/6D, EFNA1/3/4,* PSMA5/7, PSMB1/3/5, PSMC1/3/4/5, PSMD1/2/4/7/8/9/11, UBE2G2/J2, MAPK1/3,* and 21 members of ZNF family.

A substantial number of DEGs were connected with protein processing, which resulted in a significant overrepresentation of the terms "proteasome", "protein export", "arginine and proline metabolism", and "protein processing in endoplasmic reticulum" (Fig. 2A, Fig. S2A–C).

Overall, the main affected pathways can be divided into three subcategories: "cellular stress response", "metabolic adaptation", and "protein processing and degradation" (Fig. 2A).

The 164 transcripts assigned to the identified top 20 KEGG pathways were extracted from the sequencing data set, and to assess trends in up- and down-regulation of these DEGs, log_2_FC expression for pairwise comparisons of each time point was analyzed. Overall, comparing the time points revealed similarities in the HG4 vs. LG and HG4 vs. HG1 profiles, as well as their contrast to the HG1 vs. LG profile (Fig. 2B). Most (109 out of 164) transcripts were upregulated in HG1 compared to LG, whereas 133 mRNAs were downregulated in HG4 versus LG and the same amount versus HG1. Notably, most transcripts following this trend belonged to the “protein processing and degradation” and “metabolic adaptation” groups. Meanwhile, the mRNAs from the “cellular stress response” showed mixed expression patterns.

Regarding sterol biosynthesis and terpenoid backbone synthesis, the dietary changes impacted multiple genes corresponding to proteins in these pathways (Fig. S3A and B). Additionally, among the statistically significant regulated genes, several mRNAs were identified as being associated with lipid metabolism. Their expression patterns aligned with the sterol/terpenoid backbone biosynthesis genes, showing increased expression at the HG1 time point and reduced expression by the HG4 time point (Fig. 2C), including the key genes associated with these pathways ( Fig. S4).

To complement transcriptomic results, we profiled SCFAs. Corresponding to the changes in the expression of metabolic genes, in rumen fluid, acetate concentration dropped in HG4 compared to HG1 (Fig. 3A), propionate levels remained unchanged (Fig. 3B), n-butyrate concentration increased in HG1 compared to LG and HG4 (Fig. 3C), and n-valerate increased in HG1 compared to LG (Fig. 3D).Fig. 3. Metabolic changes associated with rumen adaptation. A–M Ruminal short-chain fatty acid (SCFA; A–D) as well as other indicated metabolite concentrations (E–M) across the selected time points. N Biplot showing the first canonical variate for CCA analysis of merged metabolic and transcriptomic data. O–P Gene-metabolite correlation network (|r| > 0.7, FDR < 0.05) for time point HG1 (O) and HG4 (P). R Red lines indicate negative and green lines positive correlation. The results of gene ontology “biological process” analysis for genes correlating with valerate (|r| > 0.7 and FDR < 0.05)

Expecting the most pronounced metabolome changes to occur shortly after the dietary switch between LG and HG1, we performed additional metabolomics analyses at these time points, including assessments of carboxylic acids, amino acids, biogenic amines, bile acids, and nucleotides. The concentration of detected compounds is summarized in Table S4. Notably, significant variations were found in amino acids and their derivatives. The observed changes in branched-chain amino acids (leucine and isoleucine) (Fig. 3E) and phenylalanine (Fig. 3F) might suggest altered protein turnover, aligning with the proteasome and protein degradation pathways identified in transcriptomics profiles. Additionally, an increase in taurine concentration (Fig. 3G), which is not a protein-building amino acid but is closely connected with bile acids, was detected alongside elevated levels of the bile acid deoxycholic acid (Fig. 3H), corresponding to changes in steroid metabolism indicated by the transcriptomics data.

Cadaverine, a lysine-derived diamine of microbial origin, also showed a significant decrease at HG1 (Fig. 3I), indicating an impact on microbiota-mediated amino acid metabolism. Likewise, changes in the concentration of 3-phenylpropionic acid (Fig. 3J) and phenylethylamine (Fig. 3K), both products of microbial breakdown of aromatic amino acids, suggest altered gut microbial activity and possible signaling functions.

The levels of citric acid and α-ketoglutaric acid, essential components of microbial and host TCA cycles, were altered in rumen fluid in HG1 compared to LG (Fig. 3L and M). These changes likely indicate shifts in microbial energy metabolism resulting from diet change.

Next, a correlation analysis was conducted between the 995 DEGs from the sequencing dataset and 66 metabolites. The biplot visualization of the CCA analysis revealed that most gene-metabolite pairs are located close to the origin, while a subset of variables strongly influenced the correlation (Fig. 3N). Although the canonical loadings were relatively modest, they identified a set of genes and metabolites that contribute most to the shared variation (marked in the figure), likely reflecting biologically relevant relationships. Additionally, a scatterplot of sample scores for the first canonical variate from RNA sequencing and metabolite datasets showed a very strong linear correlation, with r = 0.987 (Fig. S5A). These combined weighted patterns create a robust, coherent signal across samples, implying a strong biological link between transcriptomic and metabolomic profiles.

The visualization of the strongest correlation relationships showed a negative link between the SERTAD2 gene and several amino acids, a positive link between SLFN11 and five metabolites, and a connection between n-valerate and multiple genes primarily involved in metabolism (Fig. 3O, Table S5). Additionally, most of the other identified connections were between individual metabolites and genes.

Since changes in metabolite composition in the rumen fluid can immediately affect the rumen epithelium but also result in long-term adaptations, we aimed to assess if metabolite shifts at early time points correspond to gene expression at later time points. Among the correlation patterns with |r| > 0.7 and FDR < 0.05, most involved i-valerate or n-valerate (Fig. 3P, Table S5). These two metabolites showed positive correlations with five genes (COL1A2, COL3A1, GLT8D2, TGFB2, and SPARC). Upon submitting the list of valerate-correlating genes to GO Biological Process analysis, the main processes identified related to tissue rebuilding, indicating long-term effects (Fig. 3R).

Transcriptomics and metabolic adaptations were accompanied by morphological changes (Fig. 4A). We measured increased rumen wall thickness (RWT, Fig. 4A) as well as epithelium (Fig. 4B) and stratum corneum thickness in the papillae (Fig. 4C) in HG4. However, other histological parameters, including villi width, papillae width and length, stratum basale and stratum spinosum with stratum granulosum (Fig. S5B–F) were not affected. Corresponding to the observed changes, the expression of the IGFBP3 gene, which is associated with tissue growth, was elevated in HG4 compared to HG1 (Fig. 4E).Fig. 4. Structural changes associated with rumen adaptation. A Representative histological sections of rumen papillae at each time point (H&E stain). B–D Quantification of epithelial parameters: stratum corneum thickness, epithelial thickness, and rumen wall thickness (RWT). E Normalized count for IGFBP3 expression across time points. Data are presented as mean ± SD. ^*^P < 0.05, n = 10 or 11. Depending on the data distribution, ANOVA with Tukey’s or Kruskal–Wallis with Dunn’s post-hoc tests were applied to verify statistical significance

SARA-susceptible animals display distinct inflammatory and metabolic gene signatures persistent across diet changes

When splitting the dataset into SARA S and R animals, to characterize the group-specific phenotype, we focused on the most significant differentially expressed RNAs (P-value < 0.01) with pronounced expression changes (|log_2_FC|> 1.5) at each time point. The cut-off criteria resulted in a list of 96 DEGs.

Hierarchical clustering analysis of the DEGs demonstrated clear segregation according to SARA susceptibility (Fig. 5A). The heatmap revealed two distinct expression patterns across all three time points (LG, HG1, and HG4), with one group of genes upregulated in SARA S animals and another group upregulated in SARA R animals. This consistent pattern of differential expression, despite dietary changes, which was already evident at the LG baseline and persisted through the dietary challenge, suggests a stable, intrinsic transcriptional signature associated with SARA susceptibility.Fig. 5. Transcriptomic signatures of SARA susceptibility in rumen papillae. A Heatmap of 96 DEGs between SARA-susceptible (S) and SARA-resistant (R) cows across all time points (P < 0.01, log_2_FC > ± 1.5). B PCA trajectory plot of SARA S and R cows across three time points. Individual trajectories (thin lines) and group means (bold lines) are shown

The average PC trajectories for groups S and R show separation along PC1, indicating some differences in gene expression patterns between the two groups (Fig. 5B). This separation was observable also at the LG time point, further indicating that these transcriptional differences are not solely a consequence of the HG diet but reflect a pre-existing molecular predisposition. Statistical validation using PERMANOVA confirmed significant effects of both SARA group (R^2^ = 0.439, P < 0.001) and time point (R^2^ = 0.561, P < 0.001), with a significant group × time interaction (P < 0.001) indicating distinct temporal trajectories.

The difference along the PC1 axis between average S vs. R at LG as well as average S vs. R at HG1 is approximately 20 units, picturing a divergence in gene expression between the groups, especially considering that PC1 explains 42.62% of the variance. The separation between S and R along PC2 was less pronounced, with both groups showing more overlapping trajectories.

Moreover, both groups show changes in gene expression over time, as indicated by the movement of trajectories across the three time points. Analysis within each group revealed strong time effects in SARA R (R^2^ = 0.509, P < 0.001) and moderate effects in SARA S (R^2^ = 0.305, P = 0.046), supporting the visual observation of temporal dynamics. Group S exhibited moderate variability, with most subjects clustering around the average trajectory. In contrast, the individual trajectories within group R show higher variability, with some subjects diverging significantly from the average, indicating greater heterogeneity in gene expression patterns within this group.

Pathway analysis of the DEGs between SARA S and SARA R animals revealed three main functional categories, each with several subcategories (Fig. 6). The category "inflammation" included genes involved in immune response and inflammatory processes. The second category, "cell cycle & cell structure", was further divided into three subcategories: cell cycle, cell structure, and cell signaling & communication. The third major category, "metabolism", encompassed genes related to nutrient transport, protein metabolism, and lipid metabolism.Fig. 6. Functional categorization of genes differing between SARA phenotypes. Bar graph summarizing pathway enrichment analysis of DEGs between the SARA S and SARA R groups, classified into inflammation, metabolism (nutrient transport, lipid, and protein metabolism), and cell cycle/cell structure. The shape of the figure indicates the time point, and its size reflects the P-value for each time point

The majority of inflammation-related genes were differentially expressed between the S and R groups at the LG and HG1 time points (Fig. 6 and Fig. S5G). Similarly, protein metabolism was affected at these time points. While nutrient transport was the only term showing more statistically significant regulations at the HG4 time point.

Of particular note, CCDC196 was one of two transcripts that were statistically significantly regulated across all time points between the R and S groups. Moreover, it was the gene with the highest fold change (log_2_FC = −8) of all SARA-affected RNAs. For this reason, despite its yet unclear role, CCDC196 was included in the data visualization (Fig. 6).

The second gene to be regulated across all time points was MYO7B, with a log_2_FC ranging from −3.8 to −4.7. The gene codes myosin VIIB, which hints at differences in cell structure between the rumen of S and R groups (Fig. 6).

Network analysis of the inflammatory genes differentially expressed between SARA Sand SARA R animals, using STRING, revealed multiple documented interactions (Fig. 7A) and demonstrated significant interconnectivity, suggesting the coordinated regulation of inflammatory responses in the rumen epithelium. The network included key regulatory molecules, cytokine signaling components, and inflammatory mediators, including four members of the complement system (C7, CA3, CFH, and C5AR1).Fig. 7. Inflammatory and metabolic differences between SARA phenotypes. A STRING protein–protein interaction network for DEGs related to inflammation between SARA groups. B SCFA concentrations in rumen fluid from SARA S and R animals at the indicated time point. C–F The indicated metabolite concentrations at LG and HG1 time points. Data are expressed as mean ± SD; ^*^P < 0.05

To understand the metabolism-related gene expression changes, metabolomics data were analyzed considering SARA susceptibility. The transcriptomics changes were accompanied by differences in total SCFA, acetate, propionate, and n-valerate levels in the rumen fluid between S and R groups (Fig. 7B), while n-butyrate concentration did not differ between the groups. The majority of the SCFA-related observed differences occurred at the time points between HG1 and HG4.

In addition to SCFA, cadaverine and ribose were also differentially regulated between SARA groups (Fig. 7C, D). Additionally, upon transition from LG to HG1, deoxycholic acid and citrate levels changed in the R group, but remained unchanged in the S group, suggesting a divergence in metabolic adaptation (Fig. 7E, F).

Next, the correlation between transcriptomics and metabolomics datasets was analyzed for SARA groups. The Spearman correlation coefficients (Δρ, Fisher Z-transformed) illustrated differences between the two groups with asymmetry in the volcano plot (Fig. 8A). On the left, Δρ values drop below −2.5, some approaching −3, while on the right, few points exceed + 2, indicating stronger transcript-metabolite correlations in group R compared to group S, both in magnitude and statistical significance.Fig. 8. Correlation of metabolome and transcriptome in SARA groups. A Volcano plot indicating Spearman Rho and P-value for correlation of genes and metabolites for R and S groups. B Correlations specific to R are on the left and to S on the right side of the graph. Chord diagram visualizing gene-metabolite correlation with |Δρ|> 0.5 and P < 0.05 threshold. C Red lines indicate connections for group S and green for group R. Pathway enrichment analysis using WIKI pathways comparing DEGs correlating with metabolites (|Δρ| > 0.5 and P < 0.05) in R and S groups

The top 100 correlations with a P-value < 0.05 were visualized using a Chord diagram (Fig. 8B). These correlations showed that multiple metabolites—including amino acids (tyrosine, arginine, lysine, serine), nucleosides and related compounds (cytosine, adenosine, ribose), TCA cycle intermediate (citrate), SCFAs and microbial metabolites (butyrate, i-butyrate, n-valerate, 3-phenylpropionic acid), and fatty acid metabolism-related propionylcarnitine—as well as DEGs, were predominantly associated with the R group. Notable R-specific DEGs with multiple correlations included SEH1L, CNOT9, TMX1, and PDXDC1. Fewer DEGs, such as KCP and DNASE1, showed exclusive correlations with the S group.

All DEGs with a correlation |Δρ| > 0.5 and P < 0.05 were subjected to analysis using WIKI pathways. Although several affected pathways were shared between the S and R groups, the number of genes involved and the P-values differed, especially for proteasome degradation and cholesterol biosynthesis (Fig. 8C). This suggests a shift in tissue response to metabolite changes in SARA-susceptible animals.

Discussion

Biphasic transcriptional, structural, and metabolic adaptation of rumen papillae to an HG diet

The presented study reveals distinct temporal phases in the adaptation of rumen papillae to the HG diet. The initial response at the HG1 time point is characterized by limited transcriptional adjustments. In contrast, the pronounced gene expression changes observed by four weeks (HG4) represent a more comprehensive adaptive response, suggesting a biphasic adaptation process, aligning with previous studies [10]. The initial response typically involves physiological adjustments to maintain pH homeostasis, while longer-term adaptation includes structural remodeling of the epithelium, including thickening of the rumen wall, as shown in this study, and increasing papillae length to optimize the absorption surface [25].

Our data confirm and extend previous findings on the involvement of cholesterol and lipid homeostasis genes in rumen epithelial adaptation to the HG diet [5]. Besides validating previously identified genes, including 3-hydroxy-3-methyl-glutaryl-coenzyme A reductase (HMGCR), acetyl-CoA acetyltransferase (ACAT2), isopentenyl-diphosphate delta isomerase (IDI), and farnesyl-diphosphate farnesyltransferase 1 (FDFT1) [5], our analysis further expands the list of affected cholesterol and lipid metabolism-related factors and additionally shows the involvement of the terpenoid backbone synthesis pathway.

In a report by Steele et al*.* [5], sterol-related gene expression initially increased during the first week, followed by a decrease, suggesting early membrane remodeling with subsequent stabilization after the rapid growth phase is complete. These metabolic changes correspond to the observed structural adaptations in the stratum corneum and rumen wall thickness, which were previously reported under acidosis [26]. Alike, protein processing and degradation pathways—strongly represented in our dataset—follow this same temporal pattern, suggesting a sequence of protein turnover followed by stabilization of a new proteome composition, similar to adaptive responses observed in other tissues [13, 27]. It is also important to consider that the changes in the concentration of leucine, isoleucine, and phenylalanine, as well as microbiota-derived amino acid metabolites (cadaverine, 3-phenylpropionic acid and phenylethylamine) in the rumen fluid, impact the availability of amino acids in the rumen epithelium, potentially contributing to the protein-related gene expression changes.

Changes in ruminal SCFA concentrations appear to drive and coordinate many of the adaptive responses. We observed increased n-butyrate and n-valerate concentrations in HG1, followed by a substantial reduction in n-butyrate during HG4. Notably, butyrate serves as a preferred energy source for rumen epithelial cells and regulates numerous cellular processes, including proliferation, differentiation, and apoptosis, thus affecting barrier function and papillae growth [28–30]. One of the mechanisms via which butyrate stimulates tissue growth is downregulating IGFBP3 [31], which aligns with the mRNA level changes we measured.

Finally, we observed that some metabolites in HG1 correlate with gene expression in HG4, suggesting that metabolites may drive long-term changes and play a role in epithelial adaptation. The correlation was especially notable between valerate and a set of genes involved in tissue remodelling. Among others, these genes control extracellular matrix reconstruction (FBN2,* ELN* [32], SPARC [33],* ASPN* [34],* ADAM12* [35],* ADAMTS9* [36]), collagen deposition (COL1A2,* COL3A1*, and* COL5A3*), and TGF-β signalling, all of which are essential during epithelial restructuring. On the molecular level, valerate as a natural histone deacetylase (HDAC) inhibitor [37] may enhance the expression of the TGF-β pathway, which may promote controlled collagen deposition and tissue remodeling [38], a process essential for rumen wall thickening during HG diet adaptation. Simultaneously, valerate's activation of G-protein coupled receptors FFAR2/FFAR3 may fine-tune inflammatory responses through MAPK signaling [39], creating an environment where structural reinforcement occurs without excessive inflammation. Therefore, valerate may act as a signalling molecule that contributes to the structural reprogramming of the rumen epithelium under long-term dietary pressure. Notably, since the levels of valerate differ between SARA groups, it may impact the adaptation to the diet change.

This coordinated response involving sterol biosynthesis, protein processing, SCFA metabolism, and structural adaptation illustrates the complex, temporally regulated process by which the rumen epithelium adapts to dietary changes, beginning with initial remodeling and followed by the stabilization of a new homeostatic state.

Cellular stress response in dietary adaptation

The enrichment of glutathione metabolism, MAPK signaling, and neurodegenerative disease-associated pathways implies the involvement of cellular stress in the adaptation process. The pathways share fundamental mechanisms associated with proteostasis, oxidative stress management, and cell survival signaling. The glutathione metabolism pathway, in particular, suggests the activation of antioxidant defenses, likely in response to increased fermentation acids and the production of reactive oxygen species under HG conditions [6]. While the enrichment of MAPK signaling indicates the activation of stress-responsive cellular communication networks that coordinate adaptive responses.

These stress signatures are transient, reflecting effective cellular adaptation instead of lasting damage. This pattern aligns with the concept of the “adaptive stress response”, where the initial activation of stress pathways triggers longer-term adaptive mechanisms [40].

Despite the involvement of cellular stress and disease-related pathways, a very limited number of genes directly related to inflammation were affected by HG; thus, no inflammation-related terms were detected in our analysis. A prior study by Dionissopoulos et al*.* [41] concluded that switching to an HG diet does not affect inflammatory gene expression. The study utilized qRT-PCR and thus focused on a small selection of target genes, none of which corresponded to the genes we found regulated.

The identification of KEGG pathways such as ‘Pathways of neurodegeneration—multiple diseases’ and ‘Herpes simplex virus 1 infection’ is a common bioinformatics finding in which shared molecular machinery across biological processes leads to overlapping gene signatures. While these pathways are annotated based on neuronal or viral contexts, the DEGs driving these assignments are not neuron-specific but are instead universal regulators of core cellular processes, all of which are expressed within epithelial tissues. These include protein ubiquitination and proteasomal degradation (PSMA5/7,* PSMC1/3/4/5*,* PSMD1/2/4/7/8/9/11*,* UBE2G2*), which control protein turnover and are critical for epithelial cell cycle progression and stress response [42]. Similarly, the list included transcription regulators (e.g., multiple ZNF genes), cytoskeletal components (DNAH2/12, TUBA1D), and signaling molecules (MAPK1/3, RAC1, TNFRSF1A) that govern epithelial proliferation, migration, and differentiation [43, 44]. Their coordinated expression in the rumen epithelium likely indicates a broad stress response and tissue remodeling in response to dietary change, involving protein turnover, metabolic adjustments, and active immune responses [45]. Therefore, we interpret this result as a part of a broad epithelial stress response and tissue remodeling program triggered by the dietary change, rather than an indication of neuronal involvement.

Transcriptional and metabolic markers of SARA susceptibility

Our findings reveal a SARA-susceptibility-related transcriptional signature in rumen papillae that was persistent across all time points, including the baseline period. This indicates that the molecular predisposition to SARA is inherent and exists prior to the HG challenge. Differences in cell cycle and structure genes indicate that SARA R animals may possess a greater capacity for epithelial remodeling and for maintaining barrier integrity during dietary transitions [46]. Moreover, the expression of intercellular adhesion, cytoskeletal dynamics, cell migration, and cell turnover genes in cattle rumen has previously been associated with feed efficiency, indicating that nutrient transport and cell cycle are likely intertwined [47]. Similarly, acidosis-positive and negative rumen samples were reported to differ by expression of filament organization and cell adhesion, together with transport and absorption-related genes [48]. Correspondingly, one of the genes consistently under-expressed in SARA S individuals, MYO7B, codes a critical motor protein that plays a fundamental role in maintaining the villi's structural integrity and organization, supporting epithelial cell microvillar structures and cellular polarity [49, 50]. Its role includes intracellular trafficking, vesicle movement, and apical membrane organization. Consequently, it directly impacts the efficiency of nutrient absorption and digestive processes in intestinal epithelial cells. Since the loss of MYO7B disrupts intermicrovillar adhesion and brush border [51] and HG diets can compromise epithelial barrier function initially before compensatory mechanisms restore integrity [52, 53], the downregulation of MYO7B expression likely contributes to SARA pathogenesis. Notably, the consistently lower expression of MYO7B over all time points in SARA S animals suggests a pre-existing impairment in epithelial structural integrity, which may reduce resilience to pH challenges. This structural defect also reduces the epithelium's capacity to sense and transmit signals from the rumen lumen, such as shifts in SCFA concentration; and it allows increased translocation of microbial pathogens and their molecular patterns (e.g., lipopolysaccharide), which can activate innate immune receptors in the underlying mucosa [4, 5]. This mechanism provides a direct link between the pre-existing structural weakness and the heightened baseline inflammatory gene expression profile we observed in SARA-susceptible animals. This aligns with findings from Steele et al. [13], who demonstrated that compromised barrier function precedes SARA development.

Moreover, impaired barrier integrity may trigger local and systemic inflammatory responses that contribute to the pathogenesis of SARA [1], and heightened baseline inflammatory responses predispose the rumen epithelium to greater injury during acidotic challenges [54]. Our observation of differentially expressed inflammatory genes between SARA S and R animals supports this concept, suggesting that predisposition to SARA may partly reflect preexisting, inherent differences in the capacity to maintain barrier function under challenging dietary conditions. Furthermore, the interconnected inflammatory network that we identified highlights potential molecular targets for predicting or modifying SARA susceptibility.

The potential disparities in microbial makeup and metabolic adaptation between SARA phenotypes were reflected in the rumen levels of acetate, propionate, n-valerate, and total SCFA. Although these variances were not yet present at HG1, they were already resolved by the HG4 time point and primarily observed two weeks post-diet transition. This indicates that subsequent investigations into SARA could gain more metabolically relevant insights by focusing on week two of HG feeding. On the transcriptional level, differences in the expression of several nutrient transporters, along with genes involved in lipid and protein metabolism among SARA groups, have been measured, likely reflecting the challenges in adapting to the diet change in SARA S animals. However, based on the SCFA profiles, comparison to other time points may reveal more pronounced differences regarding the feed efficiency.

Packett and Fordham [55] showed that the utilization of citrate in rumen fluid is exclusively microbial, not host-driven, and that the addition of citrate increases SCFA production. Thus, in our experiment, the reduction in citrate in the rumen fluid of SARA R animals likely reflects enhanced microbial uptake and conversion to SCFA, which may contribute to fermentation stability. In contrast, SARA S cows did not show this citrate decline, consistent with their lower SCFA production. This suggests impaired recruitment or activity of citrate-fermenting microbial populations. Moreover, citric acid lowers rumen-fluid pH [55], implying that impaired microbiota-mediated citrate metabolism contributes to the reduction in pH in SARA S animals.

Secondary BAs, such as DCA, arise from microbial dehydroxylation of primary bile acids [56], reflecting microbial BA-metabolizing activity rather than hepatic secretion. Thus, elevated ruminal DCA in SARA R cows after HG feeding likely indicates robust microbial bile acid conversion and microbiome stability. In SARA S cows, DCA also increased, but large within-group variability impacted statistical analysis. This may reflect unstable or inefficient BA-transforming microbial populations affected by acid stress and dysbiosis, as low pH disrupts microbial communities and suppresses bile acid dehydroxylation enzymes [57]. Inconsistent DCA responses suggest impaired adaptive capacity of microbial BA pathways in susceptible animals.

Furthermore, secondary BAs modulate microbial habitat via antimicrobial and signaling functions, influencing lipid and amino acid metabolism through FXR and TGR5 pathways and affecting fermentation and host-microbe interactions [57]. Thus, DCA may represent a metabolically informative biomarker of rumen microbial adaptation to HG diet, with stable increases in resistant cows supporting efficient microbial bile acid conversion in SARA resilience.

Moreover, the correlation of the metabolome with RNA sequencing revealed additional pattern mismatch in response to the diet change between the SARA groups. The SARA R group showed a higher number of correlating gene-metabolite pairs, implying that the regulation or integration between transcripts and metabolites may be more closely coordinated or systematically altered in group R, possibly reflecting unique biological processes or network dynamics in this group. Consequently, SARA R group may experience an improved synchronization of the epithelial response to the diet- and microbiome-related metabolome, potentially facilitating better adaptation. In contrast, the S group displayed weaker and more isolated correlations, potentially indicative of dysregulated or less flexible transcriptional responses.

Finally, our analysis revealed that the CCDC196 gene shows the highest fold change (log_2_FC = −8 in SARA S vs. R group) expression difference between SARA susceptible and resistant animals. The gene was consistently affected across all analyzed time points, indicating that this is a stable feature of the susceptibility phenotype rather than an acute response to the diet. The CCDC196 gene encodes a protein of unknown function; therefore, very little is known about it. One of the very few studies mentioning it indicated that the CCDC196 gene showed the largest fold change between high- and low-fertility cows [58]. Additionally, the CCDC196 gene has been reported to be upregulated in endometrial samples of beef heifers after mating [59]. However, our report first finds a connection between this gene, the rumen, and SARA.

Together, these mechanisms likely create a compound susceptibility whereby structural, metabolic, inflammatory, and neuronal differences, present before the dietary challenge, collectively contribute to reduced adaptation capacity when facing HG diets, providing insight into potential mechanisms underlying SARA.

Study limitations and research directions

While our study provides comprehensive molecular insights, several limitations should be acknowledged. The analysis focused on three time points, potentially missing transient responses occurring between these recordings. Additionally, the functional consequences of the observed gene expression changes require validation through protein and metabolite analyses.

Further exploration of the identified stress response pathways may reveal novel strategies for mitigating the negative impacts of HG diets on rumen health. Finally, integrating microbiome data would provide a more comprehensive understanding of host-microbe interactions during dietary adaptation.

Conclusions and potential applications

In conclusion, our findings reveal an adaptive response involving the coordination of several processes culminating in metabolic adaptation and epithelial remodeling. The temporal expression pattern is characterized by biphasic transcriptional responses, including an initial stress-protective phase followed by the establishment of a new functional equilibrium that is adjusted to the HG environment. These findings enhance our understanding of rumen epithelial biology and may inform strategies to optimize rumen adaptation during dietary transitions in production settings.

Identifying transcriptional signatures associated with SARA susceptibility offers potential for developing biomarkers to identify at-risk animals, with CCDC196 and MYO7B representing promising candidates. Early identification of SARA-susceptible animals could enable targeted management strategies, including individualized feeding approaches, to mitigate the risk of acidosis [60].

Additionally, the molecular pathways identified in our study provide potential targets for nutritional interventions aimed at enhancing adaptation and reducing the risk of SARA. Supplements that support epithelial barrier function, modulate inflammatory responses, or enhance cellular stress resistance could be developed based on the pathways identified here.

Supplementary Information

Additional file 1: Fig. S1. Quality control and sample distribution for gene expression data. Fig. S2. Protein-related KEGG pathways. Fig. S3. Steroid-related KEGG pathways. Fig. S4. Expression of key genes involved in steroid and lipid metabolism. Fig. S5. Expanded structural, metabolic, and SARA-related comparisonsAdditional file 2: Table S1. Diet composition.Additional file 3: Table S2. Comparison of gene expression in rumen papillae at LG, HG1, and HG4 time points.Additional file 4: Table S3. Results of pathway analysis.Additional file 5: Table S4. Metabolomics. Concentrations (µg/mL) of detected compounds in the metabolomics analysis across samples. Compounds highlighted in red indicate low measurement quality.Additional file 6: Table S5. Correlation between genes and metabolites.