The predisposition of cows for a higher or lower MU concentration (HMUg and LMUg) was determined by deregressed proofs of an estimated breeding value for MU (EBVMU). The EBVMUs were estimated by vit Verden (Vereinigte Informationssysteme Tierhaltung, Verden, Germany) applying the established model for somatic cell score in milk [Vereinigte Informationssysteme Tierhaltung (VIT)., (2020, access 20.01.2020)]. The estimation was based on a milk record dataset of around 8 million lactation curves of milk urea values from Holstein cows in their first to third parities. The deregression of EBVMUs is a pre step for genomic breeding value estimation. The HMUg was determined as predisposition by positive values, LMUg by negative values. Based on this categorization, 10 HMUg and 10 LMUg cows were housed in pairs, comprising one HMUg and one LMUg cow, respectiveley. The phenotypes of MU concentration were obtained from the recent five milk records depicting HMUg = 278.1 ± 15.8 (mean ± SD) mg/L and LMUg = 181.54 ± 15.08 (mean ± SD) mg/L with comparable milk yields between the groups.

These cow groups were fed an isocaloric diet [10.1 ± 0.2 MJ metabolisable energy/kg dry matter (DM)], with either recommended (normal) CP content (NP with 157 ± 2.68 (mean ± SD) g CP/kg DM) or low CP content [LP with 139 ± 8.37 (mean ± SD) g CP/kg DM], resulting in a 4 ×5 balanced design (NPxHMUg, NPxLMUg, LPxHMUg, LPxLMUg; each n = 5). The experiment was conducted in 10 trials, with each trial comprising one cow pair (HMUg and LMUg), which was fed either the NP or the LP diet. The LP and NP rations were offered twice-daily for ad libitum intake as a total mixed ration. A detailed description of the diet, the feeding scheme as well as phenotypic analyses of N-metabolites (i.e., ammonia concentration in the rumen fluids, plasma urea, and glutamine concentration) is documented in Müller et al. (2021). After 2 weeks, the cows were slaughtered 4 h after morning milking and feeding by exsanguination after captive bolt stunning. Samples of rumen villi from the saccus dorsalis containing all rumen villus cell layers were collected, i.e., the basal membrane cells, which connect villi to blood stream and muscle tissue. The samples were rinsed, snap frozen in liquid nitrogen and stored at −80°C until RNA isolation. In addition, the rumen fluids (50 ml) were collected as a homogeneous mixture of the entire rumen content and stored at −20° until isolation of microbial DNA.

The rumen fluid samples (50 ml) were mixed and microbial DNA was extracted with the PowerLyzer PowerSoil DNA isolation kit (QIAGEN, Hilden, Germany) following manufacture's recommendations, but utilizing 800 μl rumen fluid sample (instead of 250 μl) and two additional heat incubation steps (step 1: 70°C, 10 min; step 2: 95°C, 10 min) before the bead beating procedure. The quantity of the extracted DNA was determined on NanoDrop ND-2000 (Thermo Fisher Scientific, Dreieich, Germany) and afterwards diluted to 10 ng/μl per sample. The diluted DNA extracts were utilized to amplify the V4 hyper-variable region of the 16S rRNA gene. The specific primer enclosed the V4-specific sequence, the sequencing flow cell adapter and a specific barcode for each primer as previously described (Kozich et al., 2013). The amplification was conducted in duplicate by polymerase chain reaction (PCR) with initial denaturation at 95°C for 2 min, followed by 35 cycles at 95°C for 30 s, 50°C for 60 s, and 72°C for 90 s, and a final extension for 10 min at 72°C using the 5 Prime HotMasterMix (5 Prime, Hamburg, Germany). The PCR products (16S amplicons) were checked on agarose gel and subsequently purified, normalized and pooled for sequencing using the SequalPrep normalization plate kit (Thermo Fisher Scientific, Dreieich, Germany) The sequencing of the amplicons was conducted on HiSeq2500 (Illumina, San Diego, CA) with 250 bp paired-end reads.

The sequence data were then trimmed for the adapter sequences and filtered with mothur software for sequences that exceeded 275 bp in length, which contained ambiguous base calls or excessively long homopolymers (>8 bases), as well as putative chimeras were removed [version 1.44.1 (Schloss et al., 2009)].The sequences were further globally aligned to the Silva reference database (release 138; https://www.arb-silva.de/, access 10.6.20) and subsequently aggregated into operational taxonomic units (OTU) considering a sequence identity of ≥97%. Taxonomic annotations for OTUs were obtained from the Silva database (release 138). Considering stratification by dispersion, the count data was rarefied to the sample with the lowest read depth (745,235 reads) as suggested by Mcmurdie and Holmes (2014).

Initially, very low abundant OTUs with <10 counts in all 20 cows were pre-discarded from further analyses. The microbial community in rumen fluids was analyzed with inverse Simpson diversity index (alpha diversity) at OTU level utilizing the R packages “vegan” (Oksanen et al., 2013) and “agricolae” (De Mendiburu Delgado, 2009). Means of Simpson indices of the predisposition groups (HMUg and LMUg) and the diet groups (NP and LP in HMUg and LMUg cows, respectively) were tested for significant differences by applying Student's t-test. The subsequent analyses were conducted at genus level, excluding genera with <10 counts in 5 or more cows to avoid potential noise by barely observed reads. The remaining dataset was tested for significantly differentially abundant genera (DAG) between the groups (HMUg–LMUg; NPxHMUg–LPxHMUg; NPxLMUg–LPxLMUg) utilizing Wald test in the DESeq2 R package (Love et al., 2014). For the contrast of predisposition (HMUg–LMUg) the trial was implemented as further fixed effect in the statistical model, accounting for different CP levels in the feed ration. For the differences between the diets, the contrasts within each predisposition group were analyzed and expressed as fold changes obtained from DESeq2 delogarithmized log2fold changes. Genera with p-values below 0.05 were considered as DAG.

The total RNA from homogenized rumen villi was isolated using TRI reagent (Sigma–Aldrich, Taufkirchen, Germany), followed by DNaseI treatment (Roche, Mannheim, Germany) and a subsequent column-based purification step utilizing the NucleoSpin RNA Kit (Macherey-Nagel, Düren, Germany). The RNA quantity was measured on NanoDrop ND-2000 (Thermo Fisher Scientific) and further assessed for integrity by application on agarose gel and via 2100 Bioanalyzer system (Agilent, Santa Clara, CA). The RNA library preparation was performed using the TruSeq Stranded mRNA kit according to the manufacture's recommendation (Illumina). The resulting libraries were quality checked on the 2100 Bioanalyzer system and quantified using the Invitrogen Qubit dsDNA HS kit (Thermo Fisher Scientific). Libraries were paired-end sequenced for 2 ×71 cycles on HiSeq 2500 (Illumina). The raw sequencing reads (fastq) were quality-checked utilizing FastQC, version 0.11.8 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The data preprocessing was conducted with Trim Galore v.0.6.5 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The low-quality reads (mean Q-score <20) and short reads (<20 bp) were filtered out. Adapter-like sequences at the 3'-end of sequence reads were trimmed. The remaining high quality paired-end reads comprised a read depth of 20,324,785 ± 3,234,163 (mean ± SD) per sample and were further aligned to the current reference genome, Bos_taurus.ARS-UCD1.2 (access 01. 08.2020) with an average alignment rate of 97.86% ± 0.01% (SD) using Hisat2, version 2.2.0 (Kim et al., 2015, 2019; Pertea et al., 2016). The number of reads uniquely mapped to each gene were extracted from the HISAT2 mapping results using HTSeq, version 0.12.4 (Anders et al., 2015).

The RNA dataset was analyzed for significantly differentially expressed genes (DEG) between groups—corresponding to the 16S microbiota sequence data (HMUg–LMUg; NPxHMUg–LPxHMUg; and NPxLMUg–LPxLMUg)—by utilizing Wald test in DESeq2 package (Love et al., 2014). Genes with <10 counts in 10 or more cows were excluded from analysis, leading to a dataset of 12,627 genes. The trial was considered as further fixed effect for HMUg–LMUg contrast. Genes with a Benjamini–Hochberg adjusted p-values <0.05 were considered as DEG. The DEG identifiers were extracted to investigate the biological function of the respective gene using the databases GeneCards (http://www.genecards.org), Ensembl (http://www.ensembl.org), and the published literature. Furthermore, genes revealing p-values below 0.01 (corresponding adjusted p-values <0.6) were imported into Ingenuity Pathway Analysis (IPA; Ingenuity® Systems, http://www.ingenuity.com) for enrichment analysis between the groups (HMUg–LMUg; NPxHMUg–LPxHMUg; NPxLMUg–LPxLMUg). Pathways with an adjusted p-value <0.05 were defined as significantly enriched between the groups. Significant pathway inhibition or activation was defined by Z-scores (significant inhibition: Z-score < -2; significant activation: Z-score >2) in IPA. The Z-scores between −2 and 0 or 0 and 2 were considered indicators of inhibition or activation, respectively. The pathways were further assigned to general functions after evaluation in Targeted Explorer (IPA). The pathways of cancer and disease were excluded.

Datasets of rumen microbial genera (n = 229 microbes) and gene expression (n = 12,627 genes) were initially transformed by variance-stabilizing transformation (DESeq2 R package). Subsequently, sparse Partial Least Squares Discriminant Analysis (sPLS-DA) was applied on both datasets. The sPLS-DA is incorporated in R package “mixOmics” (version 6.6.2, available at: http://mixomics.org) (Rohart et al., 2017) and achieves feature selection and dimension reduction simultaneously in dependence on a specific trait (Chun and Keleş, 2010). In this case, sPLS-DA was utilized to select a key subset of microbial genera and a key subset of genes that discriminate the microbial community and the rumen tissue transcriptome between HMUg and LMUg cows on two components each. The discrimination by the generated subsets was visualized in two plots.

The data subsets were further analyzed for microbial abundance–gene transcript correlations by utilizing Pearson statistics. Microbe–gene pairs exhibiting p-values below 0.01 were considered as significantly correlated. The results were visualized in a correlation map using R package “pheatmap” (version 1.0.12, available at: https://cran.r-project.org/web/packages/pheatmap) (Kolde and Kolde, 2015). The correlation map identified two dominant clusters (I and II). The genes attributing to the respective clusters were subsequently analyzed for significantly enriched pathways between HMUg and LMUg predisposition in IPA (adjusted p <0.05).

The microbial composition ranked at inverse Simpson values of 87.51 ± 39.13 (mean ± SD) in HMUg cows and 80.95 ± 34.24 (mean ± SD) in LMUg cows, indicating no significant differences between the predispositions (Figure 1A). With lower dietary CP level, alpha diversity numerically decreased, but not significantly, neither in HMUg nor in LMUg cows (Figures 1B,C).

The DESeq2 analysis at genus level revealed 10 DAG between HMUg and LMUg cows, which covered relative abundances between 0.0016 and 2.3668% across all cows (Table 1). Notably, the six most abundant genera (Prevotella, Rikenellaceae_RC9_gut_group, Methanobrevibacter, Christensenellaceae_R-7_group, F082_ge, and NK4A214_group) captured together 51.6% of total rumen abundances and did not differ significantly between HMUg and LMUg cows. Within the DAG, Unclassified Ruminococcae exhibited the broadest relative abundance (2.3668%) and manifested a 1.39-fold higher presence in LMUg compared to HMUg animals, whereas Butyrivibrio was significantly more abundant in HMUg compared to LMUg cows. Moreover, Lachnospiraceae_UCG-010 showed significantly higher occurrence in HMUg compared to LMUg, while Succinivibrionaceae_UCG-002 displayed a 2.24-fold higher occurrence in LMUg than in HMUg animals. However, highest fold changes were observed for Acetobacter and Monoglobus, with Acetobacter depicting 3.25-fold more abundance in HMUg and Monoglobus being 2.33-times more prevalent in LMUg animals.

For the comparison of NP and LP diets, 18 DAG were identified in HMUg and 19 DAG in LMUg cows, including four DAG that were previously identified in HMUg–LMUg comparison (Table 2). The reduction of dietary CP content was generally accompanied by high fold changes of genera abundances. For instance, MVP-15_ge abundances decreased 133-fold in HMUg, whereas Sharpea abundances increased 24-fold in HMUg and 166-fold in LMUg with LP feeding. Incertae sedis was less abundant in both, HMUg and LMUg animals, when cows were exposed to the LP diet. Furthermore, the abundances of Shuttleworthia increased in HMUg, but decreased in LMUg cows with the LP diet. Following dietary CP reduction, HMUg animals exhibited a further increase of Acidaminococcus abundances and of three genera belonging to phylum Desulfobacterota (Desulfovibrio, Desulfobulbus, Desulfobulbaceae_unclassified) with mean relative abundances between 0.0015 and 0.0993%. In LMUg animals, the significant increases of Streptococcus and Lachnospiraceae_NK3A20_group abundances were observed on the LP compared to NP ration. On the contrary, Alphaproteobacteria_unclassified and Succinivibrionaceae_unclassified abundances decreased remarkably under conditions of LP feeding in LMUg cows at 4.05- and 13.3-times, respectively. Interestingly, Acetobacter, whose abundance was significantly lower in LMUg than in HMUg cows (Table 1), showed a further 21-fold decrease in LMUg cows upon the reduction of the dietary CP content (Table 2).

In the filtered transcriptome dataset, 28 genes were identified as differentially expressed in the rumen between HMUg and LMUg cows (Supplementary material 1, p. 1). The highest fold changes were observed for Matrix Metallopeptidase 3 (MMP3) and Transition Protein 2 (TNP2), with higher transcript abundances in HMUg and LMUg, respectively. Furthermore, 166 genes passed the relaxed significance threshold (p <0.01) for enrichment analysis (Supplementary material 1, p. 2) and revealed 13 significantly enriched pathways that contributed almost exclusively (9 out of 13) to the immune system (Table 3). Antigen Presentation Pathway revealed the most prominent adjusted p-value. In addition, Interferon Signaling and Neuroinflammation Signaling pointed as two further candidates, displaying significant inhibition in HMUg. Whereas, Ubiquitin Like Modifier 15 (ISG15) and MX Dynamin Like GTPase 1 (MX1) were included in the former, MMP3 and Prostaglandin-Endoperoxide Synthase 2 (PTGS2) contributed to the latter. Interestingly, these four genes also ranked under the top 25 genes in terms of absolute fold changes among all annotated genes in the filtered dataset (Supplementary material 1). However, Serine Protease 2 (PRSS2) and UL16 Binding Protein 17 (ULBP17) displayed the most prominent fold changes.

Considering the reduction of dietary CP in HMUg cows, no DEGs were observed, but 68 genes passed the relaxed significance threshold for pathway analysis. This list culminated in the identification of Intrinsic Prothrombin Activation Pathway and MSP-RON Signaling Pathway as significantly enriched (Table 4; Supplementary material 2, p. 1). In contrast, LMUg animals showed 1,220 DEGs and 139 significantly enriched pathways when dietary CP content was reduced (Table 4; Supplementary material 2, p. 2, 3). Whereas, Oxidative Phosphorylation and Semaphorin Neuronal Repulsive Signaling Pathway were significantly activated when LMUg cows were fed the LP ration, additional 79 pathways were identified to be significantly inhibited, mainly representing pathways of immune response and energy metabolism as well as pathways that contribute to cell cycle and cell maintenance (Table 4).

The multivariate discriminant analysis (sPLS-DA) revealed a microbial subset comprising 50 genera in the first component and 5 genera in the second component separating HMUg and LMUg cows' rumen microbial communities (Figure 2A; Supplementary material 3, p. 1). Two genera contributed to both components, giving a subset of 53 genera that distinguished HMUg and LMUg cows. Eight genera of the subset were also identified as DAG in HMUg–LMUg contrast. The transcriptome subset included 100 genes, comprising 50 genes in each component (Figure 2B; Supplementary material 3, p. 2). Twelve genes thereof were uncovered as DEG between HMUg and LMUg.

The correlation analysis between the microbial and the gene subset separating HMUg and LMUg cows detected 157 significantly correlated microbial genera–gene pairs. The two pronounced clusters I and II were generated by two groups of microbes (A and B) and two groups of genes (I and II) that correlated contrarily (Figure 3; Supplementary material 3). In detail, Cluster I depicted gene Group I (27 genes) being positively correlated to microbial Group A (24 genera) and negatively correlated to microbial Group B (29 genera). Cluster II on the other hand showed negative correlations between gene Group II (73 genes) and microbial Group A, but positive correlations between gene Group II and microbial Group B.

The subsequent pathway analysis of gene Groups I and II revealed significant enrichment of Urea Cycle, Arginine Degradation and Citrulline Metabolism for Group I, which was mainly determined by Arginase 1 (ARG1) expression (Table 5). The ARG1 expression was negatively correlated with UCG-005 abundances (r = −0.64, p <0.01; Figure 3; Supplementary material 4). Furthermore, the differential expression of Glutamine-fructose-6-phosphate transaminase 1 (GFPT1) in Group I contributed to the significant enrichment of UDP-N-acetyl-D-glucosamine biosynthesis II (Table 5). Moreover, GFPT1 was negatively correlated with UCG-010_ge (r = −0.68, p <0.01), Lachnospiraceae_UCG-006 (r = −0.56, p = 0.011), and Roseburia abundances (r = −0.56, p = 0.01) (Figure 3; Supplementary material 4). A higher expression level of GFTP1 was accompanied by increased abundances of Lachnospiraceae_UCG-010 (r = 0.52, p = 0.02), Acetobacter (r = 0.56, p = 0.01), and Shuttleworthia (r = 0.55, p = 0.01), which had been pre-identified as DAG in the comparison of HMUg and LMUg cows (Figure 3; Table 1; Supplementary material 4).

Gene Group II revealed significant enrichment of four pathways that correspond to the immune system, including B Cell Development, IL-4 Signaling and Antigen Presentation Pathway, in overlap to the analysis of ruminal DEGs (Tables 3, 5). Particularly, genes of the Major Histocompatibility Complex (BOLA-B, BOLA-DOB, and BOLA-DRA) as well as Interferon Regulatory Factor 4 (IRF4) contributed to the significant enrichment of immune pathways in Group II (Table 5). BOLA-DRA transcripts were significantly less abundant in HMUg compared to LMUg animals (DEG) and revealed positive correlations with abundances of Roseburia (r = 0.53, p = 0.02) and Lachnospiraceae_FE2018_group (r = 0.52, p = 0.02) (Figure 3; Supplementary materials 1, 4). In addition, the presence of Gamma-Glutamyltransferase 7 (GGT7) and Leukotriene C4 Synthase (LTC4S) in Group II caused the significant enrichment of Leukotriene Biosynthesis Pathway (Table 5). The transcript abundance of both genes was highly correlated with the abundance of Prevotellaceae_UCG-004 (GGT7: r = −0.67, p <0.01, LTC4S: r = −0.52, p = 0.02) (Figure 3; Supplementary material 4). Moreover, GGT7 transcripts showed further a prominent negative correlation to Lachnospiraceae_UCG-010 abundances (r = −0.05, p = 0.02), whereas LTC4S transcripts negatively correlated with the occurrence of Coxiella (r = −0.58, p <0.01).

Furthermore, UCG-010_ge genus attracted interest due to its particularly high correlations with transcript abundances of KH Domain Containing 4 (KHDC4; r = −0.72, p = 0.0003), Ribosomal RNA Processing 9 (RRP9; r = 0.71, p <0.0005), Lipase Maturation Factor 1 (LMF1; r = 0.73, p <0.0003) and DNA Polymerase Delta 4 (POLD4; r = 0.73, p = 0.0002). Furthermore, POLD4 expression showed a strikingly high negative correlation with the abundance of Shuttleworthia (r = −0.74, p = 0.00017), whereas KHDC4 transcripts showed a prominent correlation with Acetobacter occurrences (r = 0.69, p = 0.0008).

Although no differences in the overall microbial diversity have been uncovered between HMUg and LMUg cows, some DAGs were identified. These DAG might (i) influence the amount of ruminal NPN in forms of NH₃ or NH4+, (ii) determine the utilization of blood urea by ureolytic activity, or (iii) influence the permeability for the transport of NH₃, NH4+ and urea across the rumen epithelium.

Specifically, LMUg animals displayed a significantly higher abundance of Succinivibrionaceae_UCG-002, which is an ureolytic genus in the rumen. Ureolytic bacteria metabolize urea to CO₂ and NH₃. The latter can be further utilized by other ruminal genera for microbial growth, which subsequently reduces the ruminal NPN-pool (Rosendahl, 2015). Moreover, the conversion of urea into NH₃ increases the urea concentration gradient between blood and the rumen lumen. This massively alters the urea transport rate, reduces blood urea concentration and promotes N-fixation into microbial protein (Rosendahl, 2015; Jin et al., 2016). Interestingly, significantly lower blood urea concentrations and lower blood urea pool sizes were observed in cows displaying the LMU compared to the HMU phenotype in a recent study (Müller et al., 2021). Furthermore, unclassified genera of the Succinivibrionaceae family have already been suspected as being more abundant in ruminants with increased N-use efficiency (Jin et al., 2016).

LMUg animals additionally exhibited significantly higher abundances of unclassified Ruminococcaceae. Various Ruminococcaceae genera correlated with the N-metabolism in goats, including Ruminococcus_2 displaying significantly higher abundances in animals with higher N utilization efficiency (Wang et al., 2019). Moreover, the Ruminococcaceae family has been suggested to play a fundamental role in amino acid and protein metabolism in the bovine rumen, and the Ruminococcaceae NK4A214 group has been attributed a role to N-recycling in cattle (Pacífico et al., 2021). Hu et al. (2019) found a positive correlation between the Ruminococcus abundances and the gap junctions in the rumen epithelium of yaks and postulated an absorptive response of the rumen epithelium to Ruminococcus metabolites. Some Ruminococcaceae genera were also reported to have urease activity (Patra and Aschenbach, 2018), and thus it might be conceivable that unclassified Ruminococcaceae increases the urea concentration gradient across the rumen wall and support the diffusion of urea in LMUg animals.

Furthermore, HMUg animals showed higher abundances of Butyvibrio, which is one of the most relevant butyrate producers in the rumen (Meehan and Beiko, 2014; Henderson et al., 2015; Kong et al., 2020). Interestingly, two Butyvibrio species, namely, B. fibrisolvens and B. proteoclasticus, were postulated as key species involved in protein degradation in the bovine rumen (Wallace and Brammall, 1985; Attwood and Reilly, 1995; Henderson et al., 2015). Due to their high proteolytic activities, they belong to the hyper-ammonia producing bacteria species (HAB) (Hartinger et al., 2018). Also, HAB species are known to negatively influence the N-efficiency of ruminants, since they cause steep increases of NH₃, which cannot be directly utilized by the rumen microbes and thus leads to enhanced NH₃ efflux into the bloodstream (Bento et al., 2015; Hartinger et al., 2018). Based on phenotypic categorization into HMU and LMU cow groups, HMU cows were indeed shown to exhibit numerically higher ruminal NH₃ concentrations than LMU cows (Müller et al., 2021). Furthermore, higher Butyrivibrio abundances were found to be correlated with lower N-recycling efficiency in beef cattle (Alves et al., 2020).

The influence of different CP levels on the rumen microbial community has been broadly evidenced (Lapierre et al., 2005; Aguiar et al., 2014; Patra and Aschenbach, 2018; Müller et al., 2021). In general, it is known that lower CP levels in the diets of dairy cows reduce rumen fluid NH₃ and NH4+, enhance the relative blood urea reflux and faster N-assimilation and N-fixation into microbial protein.

In both HMUg and LMUg animals, the reduction of dietary CP content induced significant changes in microbial abundances, which might be caused by microbial N-use competition and the general dietary nutrient availability. In this context, Sharpea abundances increased simultaneously with the decline of Incertae_sedis abundances in both, HMUg and LMUg cows, when feed-CP content was lowered. Sharpea is known to digest short-chain carbohydrates (Dias et al., 2017; Trabi et al., 2020), while Incertae_sedis (aggregated genus of the Ethanoligenenaceae and Ruminococcaceae families) primarily converts cellulose (Suen et al., 2011; Bailoni et al., 2021). The metabolic processes and the resulting microbial growth of cellulose digesters are generally slower than the growth of microorganisms, which metabolize the short-chain sugars. However, both microorganisms need to assimilate N components for their microbial growth (Russell et al., 2009). It might be conceivable that Sharpea depressed the growth of Incertae_sedis due to comparatively faster N-assimilation, when CP was lowered. Shuttleworthia attracted further interest since its abundance increased in HMUg, whereas significant abundance decreases were identified in the LMUg group. Several studies reported feed induced abundance alterations of Shuttleworthia, which were related to the dietary starch–fiber ratio and the energy level of the diet (Plaizier et al., 2017; Kotz et al., 2020; Tun et al., 2020). Furthermore, Zhang et al. (2019) postulated significantly different abundances between high and low yielding dairy breeds and identified Shuttleworthia as a genus displaying high relative abundances, which is in common with the present study. These findings promote Shuttleworthia as highly vulnerable and strongly adaptable genus, which might be of importance in characterizing the different ruminal adaption mechanisms in LMUg and HMUg on low CP levels.

Moreover, the abundances of Acidaminococcus and three genera of Desulfobacterota phylum (Desulfovibrio, Desulfobulbus, and Desulfobulbaceae_unclassified) increased in HMUg cows with the reduction of dietary CP content. Interestingly, several species of these genera have been previously categorized as HAB (Attwood and Reilly, 1995; Eschenlauer et al., 2002; Loubinoux et al., 2002). The growth of HAB is known to be suppressed by high NH₃ and NH4+ concentrations in a self-regulating manner (Sales et al., 2000). Interestingly, phenotypic HMU animals displayed significant decreases in previously high ruminal NH₃ and NH4+ concentrations following CP reduction, whereas in LMU phenotypes only moderate ruminal NH₃ and NH4+ declines were observed (Müller et al., 2021). Since genera abundances of the Desulfobacterota phylum appeared to be significantly increased exclusively in HMUg, it is conceivable that their growth was inhibited by high ruminal NH₃ and NH4+ levels in HMUg under NP diet feeding, and that the reduction of dietary N might have had a positive effect on their growth.

In LMUg cows, dietary CP reduction led to increased abundances of Streptococcus, which is one of the most important proteolytic members in the rumen fluid (Attwood and Reilly, 1995; Ling and Armstead, 1995). Streptococcus utilizes rapidly assimilable energy sources (starch), is known to possess ureolytic activity and is highly adaptable to various N-sources, which could provide a competitive advantage under N scarcity (Jin et al., 2016; Hartinger et al., 2018; Patra and Aschenbach, 2018).

Under conditions of LP diet feeding, LMUg cows further exhibited a significant abundance decline of Acetobacter, which had already been identified as significantly less abundant in LMUg compared to HMUg cows in this study (contrast HMUg–LMUg). The further massive abundance decline in LMUg cows from NP to LP feeding suggests Acetobacter an adaptable capacity to different N levels in the rumen, which may contribute to differential N-utilization in LMUg cows.

While the metabolic activity of the microbial community substantially affects the ruminal protein-N to NPN ratio, the gene expression in the rumen villi impact the ruminal milieu (e.g., active and passive transport processes) and thus the composition of the microbial community (e.g., immune modulation) (Penner et al., 2011; Steele et al., 2016). The pathway analysis of DEGs in rumen villi revealed the enrichment of immune response pathways, which were activated in LMUg compared to HMUg cows. Interestingly, it was recently reported in humans that the intestinal mucosa was able to respond to the microbial products with the targeted expression of receptors, immune genes and signaling molecules to maintain intestinal homeostasis (Tannock and Liu, 2020). In addition, the concept of “immune selection” was described as a way of maintaining host–microbe relationships in the long term as a result of host immune activity (Peterson et al., 2007). In this context, Bessman and Sonnenberg (2016) elucidated the essential role of Antigen Presentation and the Major Histocompatibility Complex, which dictated the host T-cell reactivity to biochemical products of commensal bacteria by the expression of innate lymphoid cells. Notably, the highest Z-score with significant activation in LMUg compared to HMUg cows was found for the Antigen Presentation Pathway in this study. Furthermore, the pathways of IL-4 Signaling, B Cell Development and Glucocorticoid Receptor identified in LMUg–HMUg contrast were enhanced, predominantly by Bovine Leukocyte Antigens (BOLA), which are prominent drivers of the bovine Major Histocompatibility Complex (Ellis and Ballingall, 1999; Lewin et al., 1999). Moreover, Zhang et al. (2021) reported a simultaneous slight activation of the Antigen Presentation Pathway with a significantly higher expression of PRRS2 in the rumen tissue of calves with comparatively higher N-efficiency. In fact, PRSS2 showed a severely higher expression in LMUg compared to HMUg cows. Thus, both antigen presentation and PRRS2 might constitute key candidates distinguishing ruminal N-utilization of LMUg and HMUg. Also, PRSS2 has further been associated with feed efficiency and performance traits in cattle (Abo-Ismail et al., 2013), which might be of interest since MX Dynamin Like GTPase 1 (MX1) and ISG15 Ubiquitin Like Modifier (ISG15) were similarly found to be higher expressed in more feed-efficient heifers, which simultaneously possessed more active immune systems (Paradis et al., 2015). All of these three genes displayed remarkably higher expression levels in LMUg compared to HMUg cows in this study driving the enrichment of Interferon and Neuroinflammation Signaling.

Among the significantly DEG, Matrix Metallopeptidase 3 (MMP3) and TNP2 exhibited the broadest expression differences between cow groups. Metalloproteases were inferred to influence the destruction of tissues by remodeling the extracellular matrix (Coussens et al., 2003; Abendaño et al., 2013) and MMP3 in specific has been recognized to influence the remodeling of adipose tissue in cattle (Jeong et al., 2017). It might be conceivable that MMP3 also attributes to structure processes in the rumen tissue and thus contribute to different absorption and diffusion capacities between HMUg and LMUg epithelia, regarding ruminal NH₃, NH4+, and urea fluxes. Furthermore, TNP2 has been reported in the context of male fertility (Pasquariello, 2015), but up to now, its function in the rumen of cattle as well as its possible role in N-metabolism remains uncertain. However, given the relationship of fertility traits and immune system activities in dairy cows, the differential expression of TNP2 demands further investigation whether LMUg and HMUg cows differ with respect to fertility parameters (Chebel et al., 2004; Hansen et al., 2004; Moore et al., 2005; König et al., 2008; Weiss et al., 2009; Hurley, 2014).

Although the reduction in dietary CP content did not trigger significantly differential gene expression in rumen tissue of HMUg animals, LMUg cows exhibited severe adaptations in their transcriptome, as characterized by more than 1,000 DEG and 139 significantly enriched pathways. The pathways identified contributed almost exclusively to the immune system, the energy metabolism and the cell cycle. Considering the significantly higher plasma urea pool size during NP compared to LP feeding in cows phenotyped for LMU (Müller et al., 2021), and the general role of immune cells in the rumen wall in protecting the host against infiltration of toxic compounds, the enrichment of these pathways may indicate the necessity for the rumen epithelia to balance the N-gradient and maintain immune homeostasis (Bessman and Sonnenberg, 2016; Tannock and Liu, 2020; Müller et al., 2021). In regard to the severe energy requirements of the immune system (Gleeson et al., 2004; Kvidera et al., 2017), it might further be speculated that the observed effects on energy metabolism pathways in LMUg cows fed the LP ration result from an altered energy demand to stabilize the ruminal N balance.

The correlation analysis between microbial and host gene subsets reflecting the rumen microbiome—host interaction, which distinguished the rumen profiles of HMUg and LMUg cows, identified two main clusters. The pathway analysis of Gene Group I revealed the enrichment of Urea Cycle and Arginine Metabolism, primarily due to increased expression of Arginase 1 (ARG1) in LMUg compared to HMUg cows. Arginase is an enzyme that catalyzes the last reaction step in the urea cycle by facilitating the synthesis of arginine into ornithine and urea. In general, the urea cycle is performed predominantly by hepatocytes in the liver (Emmanuel, 1980; Razmi et al., 2005). However, Beck et al. (2009) stated the expression of ARG1 also in the rumen epithelium of cattle, which implies catabolism of N-metabolites via the urea cycle in the rumen wall. Similarly, Emmanuel (1980) identified low ARG1 expression in the rumen epithelium of sheep. Considering the conversion of NH₃ into urea, the prominent negative correlation between ARG1 transcripts and UCG-005 abundances calls for further research. Although UCG-005 has not been mentioned directly in the context of N-efficiency in ruminants yet, Amat et al. (2021) predicted UCG-005 a negative impact on Methanobrevibacter abundances in yearling heifers, which is known to contain various N-fixating strains (Poehlein et al., 2018).

Further insights into a possible differentiation of HMUg and LMUg were displayed by Gene Group II, which revealed the enrichment of four immune system pathways, particularly attributed to the expression of BOLA genes. Interestingly, BOLA polymorphisms were not only reported to determine inherent and acquired immune defenses (Rupp and Boichard, 2003), they were also capable to modulate the composition of colostrum microbiota in dairy cows (Derakhshani et al., 2018). Furthermore, BOLA gene expression is known to influence the symbiotic microbial community in the gastrointestinal tract of humans (De Palma et al., 2010), mice (Toivanen et al., 2001), and fish (Bolnick et al., 2014). In the present study, BOLA-DRA transcripts were significantly more abundant in LMUg compared to HMUg cows and showed further striking positive correlations to Roseburia abundances. The contribution of Roseburia to immune responses and to the N-metabolism in cattle, goats, and rabbits has been postulated by various authors (Wang et al., 2011; Alves et al., 2020; Sun et al., 2020). The observed higher abundances of Roseburia and BOLA-DRA transcripts in LMUg compared to HMUg cows might picture a microbial–host interplay that contributes to the distinction of ruminal patterns between HMUg and LMUg cows by affecting immune response and N-metabolism-associated processes.

In addition to the correlation with Roseburia, BOLA-DRA transcripts revealed a severe positive correlation to Lachnospiraceae_FE2018_group abundances, which belongs—like as Roseburia—to the Lachnospiraceae family. Lachnospiraceae_FE2018_group had been mentioned in terms of intestinal barrier function, antioxidant balance, and intestinal inflammation in the gut immune system of broilers (Kong et al., 2020; Liu et al., 2021). Moreover, Wang et al. (2021) postulated a significant positive correlation between Lachnospiraceae family and intestinal inflammation in dairy cows. The Lachnospiraceae family was also found to be significantly more abundant in beef steers with low N-retention efficiency and higher urinary N-excretion (Alves et al., 2020). Furthermore, high abundances of Lachnospiraceae were uncovered in the colon and caecum of goats with low N-utilizing phenotype (Wang et al., 2019). Additionally, Lachnospiraceae occurrences were found to be correlated with a low feed conversion rate and low feed efficiency in beef steers (Hernandez-Sanabria et al., 2010; Carberry et al., 2012).