Skip to main content

Microbiota-host crosstalk in the newborn and adult rumen at single-cell resolution

Abstract

Background

The rumen is the hallmark organ of ruminants, playing a vital role in their nutrition and providing products for humans. In newborn suckling ruminants milk bypasses the rumen, while in adults this first chamber of the forestomach has developed to become the principal site of microbial fermentation of plant fibers. With the advent of single-cell transcriptomics, it is now possible to study the underlying cell composition of rumen tissues and investigate how this relates the development of mutualistic symbiosis between the rumen and its epithelium-attached microbes.

Results

We constructed a comprehensive cell landscape of the rumen epithelium, based on single-cell RNA sequencing of 49,689 high-quality single cells from newborn and adult rumen tissues. Our single-cell analysis identified six immune cell subtypes and seventeen non-immune cell subtypes of the rumen. On performing cross-species analysis of orthologous genes expressed in epithelial cells of cattle rumen and the human stomach and skin, we observed that the species difference overrides any cross-species cell-type similarity. Comparing adult with newborn cattle samples, we found fewer epithelial cell subtypes and more abundant immune cells, dominated by T helper type 17 cells in the rumen tissue of adult cattle. In newborns, there were more fibroblasts and myofibroblasts, an IGFBP3+ epithelial cell subtype not seen in adults, while dendritic cells were the most prevalent immune cell subtype. Metabolism-related functions and the oxidation-reduction process were significantly upregulated in adult rumen epithelial cells. Using 16S rDNA sequencing, fluorescence in situ hybridization, and absolute quantitative real-time PCR, we found that epithelial Desulfovibrio was significantly enriched in the adult cattle. Integrating the microbiome and metabolome analysis of rumen tissues revealed a high co-occurrence probability of Desulfovibrio with pyridoxal in the adult cattle compared with newborn ones while the scRNA-seq data indicated a stronger ability of pyroxidal binding in the adult rumen epithelial cell subtypes. These findings indicate that Desulfovibrio and pyridoxal likely play important roles in maintaining redox balance in the adult rumen.

Conclusions

Our integrated multi-omics analysis provides novel insights into rumen development and function and may facilitate the future precision improvement of rumen function and milk/meat production in cattle.

Background

Ruminants, which were domesticated to provide meat and milk for humans during the Neolithic [1], largely rely on the rumen system to digest human inedible plant mass and transform this into edible protein [2,3,4,5]. Newborn (NB) ruminants are not able to chew the cud and depend on milk-based diets. They are considered functionally monogastric animals due to the underdeveloped rumen and less abundant resident microbiota [6]. In adult (AD) ruminants, the rumen displays a distinct morphological structure and has essential biological functions. For example, the rumen expands gradually after birth to make up more than 70% of total digestive tract volume for efficiently digesting forage-based diets [7, 8]. A previous study also indicated that the rumen was the organ that had the most significant transcriptional differences between newborn and adult stages by comprehensively characterizing transcription of the entire gastrointestinal tract [9]. The rumen is colonized with a complex and diverse community of microbiota. A recent study reported that the ruminal transcription of genes related to the “respiratory electron transport chain” was positively correlated with age as well as with the relative abundance of ruminal bacterial genera [10]. However, to date, there have been few studies explored the development of mutualistic symbiosis between rumen and its microbial colonizers.

The microbiota in the rumen has been classified into three major groups: microbiota in solid digesta, microbiota free-floating in the liquid fraction, and the mucosal microbiota. The microbiota in solid digesta and free-floating in the liquid fraction has been the focal point of the rumen physiology research, whereas the mucosal microbiota has received less attention [11]. Previous studies suggested that the analysis of rumen tissue-attached bacteria may provide a better understanding of ruminant metabolism and host-microbe interactions [11, 12]. To date, however, the underlying interactions between rumen mucosal microbiota and host cells remain largely unknown. It is known that almost the entire rumen epithelium is covered by microbes [13, 14] and that colonization by ruminal epithelial bacteria is age-related [15]. In this respect, the rumen is also a useful model of how microbe-host interactions develop.

Nutrient absorption and microbial interaction are facilitated by the large surface area of the rumen, which is lined with abundant papillae of keratinized stratified squamous epithelium composed of four different cellular strata: basale, spinosum, granulosum, and corneum [16]. Recent advances in single-cell RNA sequencing (scRNA-seq) have paved the way for a more granular analysis of cell types, and the finding that at least 15 cell types can be detected in the intestinal columnar epithelium of humans [17] led to speculation that more discrete epithelial cell populations exist in the rumen [18]. In this study, we have taken a multi-omic approach that includes scRNA-seq to explore the cell-type composition, associated microbiome, and metabolome of rumen epithelium in a comparative analysis of undeveloped rumen tissues from newborn calves and mature rumen tissues from adult lactating cattle. Our scRNA-seq data from the rumens of these adult cattle has also been reported and analyzed in two other studies (embarked upon after this one, but already published [19, 20]).

In one study [19], we compared scRNA-seq data and metabolites from multiple tissue types (rumen, reticulum, omasum, abomasum, ileum, rectum, liver, salivary gland, mammary gland, and peripheral blood), originating from the same lactating adult dairy cattle. In the other study [20], we investigated the contribution of interactions between adult rumen epithelial cell subtypes and rumen digesta microbiota to fiber utilization by integrating the data of meta-genome assembled microbial genomes and scRNA-seq [20]. In the current study, we used scRNA-seq for analyzing cell-type/subtype composition and function of NB and AD rumen in dairy cattle and performed the cross-species analysis between the cattle rumen and the human stomach and skin. 16S rRNA gene sequencing was conducted to analyze the mucosal bacterial composition. To gain insights into how the mucosal bacteria interact with the host’s cells, we performed metabolomics of NB and AD rumen tissues. Integration of epithelial microbiome, metabolome, and single-cell transcriptome was further explored to infer the microbiota-rumen cell-type interactions.

Results

Single-cell atlas of the rumen in NB and AD cattle

To delineate the cell-type composition dynamics of the rumen between NB and AD stages, we generated scRNA-seq profiles for rumen tissues from three newborn calves and three adult dairy cattle. After quality filtering, we profiled a total of 49,689 high-quality single cells and classified the cells into coherent transcriptional clusters (Fig. 1A). Next, we annotated the cell clusters by the average expression of canonical marker genes and identified four major cell types: immune cells (high PTPRC expression), fibroblast (high COL3A1 and SPARC expression) [21], endothelial cells (high PECAM1 and CLDN5 expression) [22], and epithelial cells (high EPCAM expression) [23] (Fig. 1A, B). Multi-donor analysis of rumen tissues showed limited effects on cell-type discovery in NB or AD tissues (Fig. 1C). Epithelial cells were the most abundant cell type in both NB and AD groups, accounting for 87.52 and 81.10% of total cells, respectively (Fig. 1D). Immune cells were mainly identified from adult cattle, while most endothelial cells and fibroblasts were from newborn calves (Fig. 1E). The number of the four major cell types was also available in Table S1. There are many subtypes of immune cells and non-immune cells; however, little is known about them in the rumen. We next analyzed the rumen immune cells and non-immune cells separately (see “Methods”).

Fig. 1
figure 1

Dissection of the ruminal cell type composition with scRNA-seq in NB and AD dairy cattle. A UMAP plots of cells from the 6 samples profiled in this study, with each cell color coded to indicate the associated cell types. B The UMAP plots show the expression of curated genes in the cell types defined. C The UMAP maps of adult (AD) and newborn (NB) rumen single-cell data, with each cell color coded to indicate the associated samples. D Pie chart showing the relative percentages of different major cell types in newborn calves and adult cattle. E The cell-type contribution bar charts show the distribution of major cell types in NB and AD rumen samples. NB: newborn; AD: adult

Immune cell subtypes and difference between NB and AD cattle

In total, immune cells were clustered into six separate subtypes (Fig. 2A): Th17 cells, MKI67+ Th17 cells, γδ T cells, natural killer T (NKT) cells, dendritic cells (DC), and plasmacytoid DC. Th17 cells (Cluster 0) expressed the classic markers, including CD3E, CD4, IL17A, and IL17F [24] (Fig. 2B; Additional file 2: Table S2). Except for the markers of Th17 cells, Cluster 4 uniquely expressed MKI67 and was annotated as MKI67+ Th17 cells (Fig. 2B). In humans, γδ T cells usually present a CD3E+CD4-CD8A- cell surface phenotype and express the TCRδ constant region-encoding segment TRDC [25]. Cluster 1 showed high expression of CD3E and a bovine homolog of TRDC but no expression of CD4 and CD8A, which suggested that these cells were γδ T cells (Fig. 2B; Additional file 2: Table S2). Cluster 3 was defined as NKT cells due to its specific expression of CD3E, CTSW, KLRK1, NKG7, CCL5, and CD8A [26, 27] (Fig. 2B; Additional file 2: Table S2). Cluster 2 was DC for highly expressing FLT3, CTS3, and BOLA-DRA as well as the conventional type 2 dendritic cell marker genes FCER1A and CD1E, and Cluster 5 was believed to represent plasmacytoid DC that expressed SPIB, BLNK, and ITM2C [28, 29] (Fig. 2B; Additional file 2: Table S2). Our other paper, which analyzed the same scRNA-seq data form the AD rumen but without the NB rumen, also identified the Th17 cells, MKI67+ Th17 cells, γδ T cells, and DC, but no NKT cells and plasmacytoid DC [19]. However, we found the NKT cells and plasmacytoid DC by integrating scRNA-seq datasets of the NB and AD rumen in this study. The reason for the inconsistency is that when analyzing solely the scRNA-seq data of the AD rumen, the NKT cells and plasmacytoid DC were algorithmically grouped with a more populous cell type due to the small number (see below). This is a predictable artefact of the annotation scheme: entire cell clusters, rather than individual cells, were annotated in each tissue. Th17 cells represented the most prevalent immune cell type in the AD rumen, while DC was the prevalent immune cell type in NB rumen. Compared with NB calves, the average percentage of Th17 cells (67.83% vs. 11.13%) and γδ T cells (22.41% vs. 4.02%) were increased in the AD cattle, while the average percentages of DC (2.40% vs. 62.44%), plasmacytoid DC (0.06% vs. 1.88%), and NKT (5.05% vs. 18.11%) cells were decreased in the AD cattle (Fig. 2C). These results provide novel insights into the different rumen immune microenvironment between NB and AD stages.

Fig. 2
figure 2

Dissection of the immune and non-immune cell composition in NB and AD rumen tissues. A UMAP plots of immune cells, with each cell color coded to indicate the associated cell types. B The violin plots showing gene expression of markers of each immune cell types. C Bar graph showing the relative percentages of different immune cell types in newborn calves and adult cattle. D UMAP plots of non-immune cells, with each cell color coded to indicate the associated cell types. E The dot plot visualization of each cell type in rumen non-immune cell data. Dot size represents the percentage of cells within a cell type, and the color encodes the expression level. F The cell-type contribution bar charts show the distribution of cell types in NB and AD rumen samples. NB: newborn; AD: adult. NKT: natural killer T cells, DC: dendritic cells; MC1: mitotic cell 1; MC2: mitotic cell 2; BC1: basal cell 1; BC2: basal cell 2; BC3: basal cell 3; SC1: spinous cell 1; SC2: spinous cell 2; cg-like SC: channel-gap-like spinous cell; GC1: granule cell 1; GC2: granule cell 2; GC3: granule cell 3; TNFRSF6B+ EC: TNFRSF6B+ epithelial cell; IGFBP3+ EC: IGFBP3+ epithelial cell; VEC: vascular endothelial cell

Non-immune cell subtypes and cross-species comparison

In addition to immune cell subtypes, we observed 17 cell clusters of non-immune cells including 13 subtypes of rumen epithelial cells (EC), 3 fibroblast subtypes, and 1 vascular endothelial cell (VEC) subtype based on the feature plots of known cell-type markers (Figs. 2D, E; Additional file 2: Table S3). Basal cells (BC), spinous cells (SC), and granule cells (GC) are the living epithelial cell types currently characterized in the rumen [30, 31]. Three BC (high KRT14 and KRT5 expression) [32,33,34], three SC (high S100A8 or KRT6A expression) [35], and three GC (high CLDN4, CLDN1, or DLK2 expression) [30, 31, 36] cell subtypes were discovered in the rumen (Fig. 2D, E). In addition, we revealed two mitotic cell (MC) subtypes based on the expression of the marker genes of KRT14 and MKI67 [34, 37] (Fig. 2E). The mitotic cell 1 (MC1, Cluster 11) also highly expressed KRT19, CENPA, CENPF, etc., and the mitotic cell 2 (MC2, Cluster 7) also highly expressed TOP2A, UBE2C, etc. (Fig. 2E; Additional file 2: Table S3). Cluster 9 spinous cell subtype specifically highly expressed GJA1, encoding a channel gap junction protein [38]; thus, it was annotated as channel-gap-like spinous cells (cg-like SC). Interestingly, Cluster 10 and Cluster 15 commonly expressed BC and GC markers but uniquely expressed TNFRSF6B and IGFBP3, respectively. Thus, Cluster 10 and Cluster 15 were defined as TNFRSF6B+ EC and IGFBP3+ EC, respectively. In our other paper [20], we also found the MC, BC, cg-like SC, SC, and GC cell types and all of them had more subtypes than that of this study (3 vs. 2, 4 vs. 3, 3 vs. 1, 3 vs. 2, and 5 vs. 3, respectively). This was in line with expectations because our other paper [20] carried out reclustering analysis on the epithelial cells collected from the same scRNA-seq data of AD rumen [19] with higher resolution (1.6 vs. 0.8), which can identify more cell subtype clusters. It should be noted that the TNFRSF6B+ EC was observed in this study but not in our other study. Given that the GC_5 cell cluster identified in our other study uniquely highly expressed the TNFRSF6B gene but not BC marker genes [20], the TNFRSF6B+ EC cell cluster may be a combination of GC_5 and basal cells in this study. This suggests that it may be necessary to perform the reclustering analysis with higher resolution on the epithelial cells identified in this study in the future.

Diverse rumen structural cell subtypes were identified, including VEC (Cluster 8), myofibroblasts (Cluster 13), PAM+ fibroblast (Cluster 14), and MFAP5+ fibroblast (Cluster 16) (Fig. 2D, E). A list of differentially expressed genes for each cluster was also available in Table S3. In our other paper [19], the myofibroblasts, PAM+ fibroblast, and MFAP5+ fibroblast cell clusters were not identified. This is because they could not form the cell cluster and algorithmically grouped with a more populous cell cluster due to the small number (only 9, 2, and 0, for the myofibroblasts, PAM+ fibroblast, and MFAP5+ fibroblast cell clusters, respectively; Additional file 2: Table S1) when analyzing solely the scRNA-seq data of the AD rumen. The results suggest that a sufficient number of cells are needed to identify a cell type in the single-cell RNA sequencing studies.

The rumen is the hallmark organ of ruminants, which is likely to have different cell composition from human and other monogastric animal stomach or other organs. To compare the cattle rumen and the human cell landscape, we downloaded a human stomach single-cell RNA-seq dataset from a previous study [23] and performed cell clustering analysis. We identified 23 cell clusters in the human stomach and identified 8 epithelial cell clusters based on highly expressed marker genes, such as EPCAM, KRT19, KRT18, KRT8, PGA5, CKN1, CHGA (Additional file 1: Figs. S2A and B; Additional file 2: Table S4). Considering that circulating immune cells that derive from the same lineage are widely distributed within the mammalian body, the current study focused on the cross-species comparison of epithelial cells (see “Methods”). A correlation heatmap among 13 rumen and 8 human stomach epithelial cell clusters showed that the similarities of epithelial cells between the two species were very low; all the epithelial cell subtypes of the cattle rumen have a weak correlation with the cell clusters of the human stomach (area under the receiver operating characteristics (AUROC) score < 0.9) (Additional file 1: Fig. S2C).

Similar to the rumen tissue, the human skin tissue is also composed of stratified squamous epithelium [38] and shared a significant number of transcriptomic features and conserved gene expression patterns with the rumen [39]. Therefore, we also downloaded the human skin single-cell RNA-seq dataset from a previous study [37] and performed cell clustering analysis. In the human skin, we identified 17 epithelial cell clusters based on highly expressed marker genes, such as KRT14, KRT5, KRT15, KRT19, KRT18, KRT6A, KRT10, KRT1, and FLG (Additional file 1: Figs. S3A and B; Additional file 2: Table S5). By performing cross-species analysis of 13 rumen and 17 human skin epithelial cell clusters (see “Methods”), we found that the gene expression patterns of most epithelial cell clusters were not conserved (AUROC score < 0.9); only 3 epithelial cell clusters of human skin had a strong correlation with rumen cell clusters (AUROC score > 0.9) (Additional file 1: Fig. S3C). Cluster 10 in the human skin had a strong correlation with cg-like SC of the rumen (Additional file 1: Fig. S3C) and expressed some marker genes of cg-like SC, such as GJA1, FABP5, and KRTDAP (Additional file 2: Table S5). Despite this, only 18 of the top 100 marker genes were shared between the two cell clusters (Additional file 2: Tables S3 and S5), and other genes of the top 100 genes in cattle that were different from humans were enriched in the “fatty acid binding” (such as S100A8, S100A9, and FABP4 genes) and “fatty acid derivative metabolic process” (such as GSTA1, MGST3, and ATP6V1B1 genes) based on the Gene Ontology (GO) enrichment analysis, which may benefit cattle’s energy acquisition, as fatty acids, rather than glucose, are the main source of energy in ruminants. Cluster 18 in the human skin had strong correlations with the MC1 and MC2 of the rumen; and Cluster 11 in the human skin had a strong correlation with the MC2 (Additional file 1: Fig. S3C). Clusters 18 and 11 in the human skin also highly expressed some marker genes of MCs, such as MKI67, CENPF, and HMGB2 (Additional file 2: Table S5). However, only 38 and 50 of the top 100 marker genes of MC1 and MC2 were shared with that of Cluster 18 in the human skin, respectively, and only 35 of the top 100 marker genes were shared between the cattle MC2 and Cluster 11 of the human skin (Additional file 2: Table S3 and S5). Taken together, our findings suggested that almost all epithelial cell types in the rumen of cattle and the stomach and skin of humans were not similar.

Although newborn calves are considered functionally monogastric animals due to the underdeveloped rumen [6], it is interesting to note that almost all the rumen epithelial cell subtypes were presented in both NB and AD cattle, except for IGFBP3+ EC, which were only identified in NB calves (Fig. 2F). The myofibroblast and PAM+ fibroblasts were predominantly distributed in the NB calves, and the MFAP5+ fibroblasts were not found in the AD cattle in this study (Fig. 2F). To some extent, the dissociation bias inherent in all solid tissues single-cell experiments may cause spurious changes in cell numbers. However, through reanalysis of previously published transcriptome data obtained by the bulk RNA-seq experiment from newborn [10] and adult [40] cattle rumen tissues, we observed that high expression levels of the fibroblast marker genes COL3A1(9,230 vs. 636.8, counts per million (CPM)), PAM (293.3 vs. 148.8, CPM), and MFAP5 (173.8 vs. 7.6, CPM) in the newborn rumen tissues whereas very low expression levels in adult rumen tissues (Additional file 1: Fig. S1).

Functional differences of cell types between NB and AD cattle

In addition to cell-type composition heterogeneity, we further deciphered the functional differences in the specific pairs of cell types between NB and AD rumen. We first divided the dataset into NB and AD groups and compared gene expression patterns of individual cell types between groups (see “Methods”). Compared to the NB group, we identified 266 upregulated and 362 downregulated differentially expressed genes (DEGs) (|log2FC| > 0.5, adjusted p-value < 0.05), both of which were shared by at least two cell types (Fig. 3A, B; Additional file 2: Table S6). We also identified the cell-type-specific upregulated and downregulated DEGs (Fig. 3A, B; Additional file 2: Table S6). Interestingly, we note that the immune cell subtypes had a larger number of cell-type-specific upregulated DEGs than cell-type-specific downregulated DEGs whereas the epithelial cell subtypes had fewer cell-type-specific upregulated DEGs than cell-type-specific downregulated DEGs in the AD cattle.

Fig. 3
figure 3

Changes in the transcriptional profiles of different cell types between NB and AD rumen. A, B Heatmaps show the distribution of upregulated (A, blue) and downregulated (B, red). The upper part (dotted lines) represents the DEGs shared by at least two cell types, and the lower panel represents the unique DEGs of each cell type, and the number below the doted box represents the number of the unique DEGs in each cell type. C, D The representative GO terms of cell types (enriched in upregulated and downregulated DEGs, respectively; p-value < 0.05). NB: newborn; AD: adult. NKT: natural killer T cells, DC: dendritic cells; MC1: mitotic cell 1; MC2: mitotic cell 2; BC1: basal cell 1; BC2: basal cell 2; BC3: basal cell 3; SC1: spinous cell 1; SC2: spinous cell 2; cg-like SC: channel-gap-like spinous cell; GC1: granule cell 1; GC2: granule cell 2; GC3: granule cell 3; TNFRSF6B+ EC: TNFRSF6B+ epithelial cell; VEC: vascular endothelial cell

To further explore which biological processes the DEGs are involved in, we performed GO enrichment analysis of up- and downregulated DEGs. The upregulated DEGs across all the rumen epithelial cell subtypes in AD cattle were mainly enriched in GO terms of the “fatty acid derivative metabolic process,” “oxidation-reduction process,” and the metabolic processes of retinoid, antibiotic, icosanoid, diterpenoid, and cofactor (Fig. 3C), suggesting dramatic shifts in physiology and metabolism of rumen between NB and AD cattle. However, the downregulated DEGs involved GO terms were not consistent among different cell types (Fig. 3D). For instance, the downregulated DEGs in cg-like SC were specifically enriched in the “regulation of animal organ morphogenesis.” The GO term “regulation of epidermis development” was only enriched in the MC1, basal cell 1 (BC1), basal cell 2 (BC2), and spinous cell 2 (SC2). Nevertheless, GO analysis showed that downregulated genes were mainly related to epithelial cell proliferation and receptor signaling pathways across many epithelial cell subtypes. Taken together, these results indicate a series of cell-type-specific molecular changes featuring cell proliferation and development at NB stage, and oxidation-reduction and nutrient metabolism at AD stage in the rumen tissues.

Rumen mucosal microbial composition difference between NB and AD cattle

Mucosal microbiota lives in close contact with host cells of rumen tissues, which can execute coordinated functions contributing to the diverse physiological processes of the rumen [41]. To further explore how mucosal microbiota interact with host cells in the newborn and adult rumen, we performed 16S rRNA gene sequencing analysis of the same rumen tissues from NB and AD cattle. To minimize the effects of sequencing depth on alpha diversity measure, the number of sequences from each sample was rarefied to 25,349, which yielded an average Good’s coverage of over 99%. Rarefaction curves approximately trended to a plateau and sequencing coverage was saturated at 15,000 reads (Additional file 1: Fig. S4). We found higher bacterial diversity (2.32 vs. 5.85, P = 0.004) and richness (136.7 vs. 770.3, P = 0.007) in the AD rumen compared with the NB group (Fig. 4A), which is consistent with the previous studies [15, 42]. The composition of rumen bacteria at the genus level also showed clear discrimination between the two stages (Fig. 4B). Furthermore, a total of 94 significantly different (LDA > 3.5, P < 0.05) bacterial taxa were identified between NB and AD groups (Fig. 4C). The NB group was only enriched with 5 bacteria genera (Fig. 4C; Additional file 1: Fig. S5), including the Megasphaera, Enterococcus, Acidaminococcus, Streptococcus, and Lactococcus. Compared with NB, the AD group was significantly enriched in 32 bacteria genera (Fig. 4C; Additional file 1: Fig. S5), including Butyrivibrio and Prevotella, which are known to produce SCFAs [15, 43]. The Desulfovibrio genus also had a higher relative abundance (1.78% vs. 0.01%) in the AD rumen (Fig. 4C), and the Desulfovibrio was also verified by the bacterial fluorescence in situ hybridization (FISH) (Fig. 4D). We further determined the bacterial densities for the Desulfovibrio genus and other genera (including Enterococcus, Butyrivibrio, and Prevotella) being differentially abundant between NB and AD groups using absolute quantitative real-time PCR. The bacterial densities of Desulfovibrio genus and Butyrivibrio genus in AD group were significantly higher than that in NB group (P < 0.01) (Fig. 4E). The bacterial density of Prevotella genus in AD group was tended to be higher than that in NB group (P = 0.08) (Fig. 4E). The bacterial density of Enterococcus genus in NB group was significantly higher than that in AD group (P < 0.01) (Fig. 4E).

Fig. 4
figure 4

Epithelial microbiota in the rumen tissue between NB and AD dairy cattle. A The difference of α-diversity indexes of ASVs between NB and AD groups. B The epithelial bacterial composition in NB and AD groups on genera level (relative abundances > 5%). C The linear discriminant analysis effect size (LEfSe) analysis and cladogram representation show that bacterial taxa differed significantly different between NB and AD groups. D The sections were stained with DAPI (epithelium, blue) and bacterial FISH probe (green) targeting Desulfovibrio demonstrating its higher abundance in AD than NB rumen tissues. E Differences of rumen epithelial bacteria density between NB and AD groups. NB: newborn; AD: adult

Metabolome profiles of rumen epithelial tissue in NB and AD animals

In the above sections, we revealed that all the rumen epithelial cell types exerted a stronger function of oxidoreductase activity, and the Desulfovibrio had higher relative abundance at AD rumen tissues compared to NB group, which could take part in the host oxidation-reduction process. In many ecosystems, the bacteria play a significant role as manifested by their dynamic interactions with the host. The microbe-host interactions are mediated by the interplay of various molecular components that are expressed by the host and the bacteria. Therefore, we performed a metabolomics analysis in the NB and AD rumen tissues (see “Methods”). In total, 718 metabolites were detected (Fig. 5A), which belong to lipids (n = 167), amino acids and their derivatives (n = 94), organic acid and its derivatives (n = 58), nucleotide and its derivatives (n = 48), alcohol and amines (n = 29), heterocyclic compounds (n = 26), benzene and substituted derivatives (n = 22), carboxylic acids and derivatives (n = 21), coenzyme and vitamins (n = 11), bile acids (n = 9), tryptamines, choline, pigments (n = 8), hormones and hormone-related compounds (n = 6), and other unnamed metabolites (n = 211). Among the 167 lipids, 31, 28, 13, 10, 8, 5, and 5% were carnitines, lysophospholipids, oxidized lipids, fatty acyls, fatty acids, unsaturated fatty acids, and glycerides, respectively (Fig. 5A). Among the 94 amino acids and their derivatives, amino acid derivatives, amino acids, and small peptides accounted for 49, 38, and 13%, respectively (Fig. 5A). A total of 112 metabolites were significantly higher (variable importance in projection (VIP) ≥ 1, |Log2FC| ≥ 1) while 179 were significantly decreased (VIP ≥ 1, |Log2FC| ≥ 1) in the AD group when compared with the NB group (Fig. 5B; Additional file 2: Table S7).

Fig. 5
figure 5

Desulfovibrio associated metabolites and cell types in the rumen. A The classification of metabolites identified in the rumen tissue. B The differential metabolites between NB and AD groups. C The Venn diagram represents the overlaps of the upregulated metabolites and the microbiota-dependent metabolites that have high co-occurrence probabilities with Desulfovibrio. D The co-occurrence probabilities of the four metabolites. E Gene scoring analysis of rumen epithelial cells between NB and AD groups using the vitamin B6 binding gene set. p values are calculated from two-sided Wilcoxon rank-sum tests. NB: newborn; AD: adult. MC1: mitotic cell 1; MC2: mitotic cell 2; BC1: basal cell 1; BC2: basal cell 2; BC3: basal cell 3; SC1: spinous cell 1; SC2: spinous cell 2; cg-like SC: channel-gap-like spinous cell; GC1: granule cell 1; GC2: granule cell 2; GC3: granule cell 3; TNFRSF6B+ EC: TNFRSF6B+ epithelial cell

Desulfovibrio interacts with rumen epithelial cells via pyridoxal

The production of small molecules is an essential route of rumen bacteria influencing host physiology [44, 45], yet it is difficult to accurately monitor the diversity of molecules produced by mucosal bacteria. To address this gap, we used the gut microbial metabolites database [46] and the microbe–metabolite vectors (mmvec) neural networks analysis [47] to identify host-microbe interactions by calculating the co-occurrence probabilities between mucosal bacteria and the ruminal tissue metabolites (see “Methods”). We found that 142 metabolites had high co-occurrence probabilities (the inferred conditional probabilities > 1) with Desulfovibrio genera, and 38 of them are microbiota-dependent metabolites against a reference library of microbiota-dependent metabolites [46] (Additional file 2: Table S8). The Venn diagram represents the overlaps of these 38 microbiota-dependent metabolites and the upregulated metabolites, which included betaine, pyridoxal, nicotinic acid, and 4-hydroxybenzaldehyde (Fig. 5C, D). As one form of vitamin B6, pyridoxal is a coenzyme participating in many metabolic reactions, including epithelium proliferation, transportation of amino acids, and inflammatory response [48,49,50]. And pyridoxal, as an electron donor for sulfate reduction, can be oxidized by Desulfovibrio [51], contributing to reducing reactive oxygen species (ROS).

Finally, we computed molecular signature scores of the vitamin B6 binding of rumen epithelial cell types between the two groups using the gene set scoring analysis. The results showed that rumen epithelial cell types scored significantly higher (P < 0.05) for vitamin B6 binding in AD group compared with the NB group (Fig. 5E), which implies that the pyridoxal may be implicated in the interaction between the Desulfovibrio and the rumen cell types. Further study to unveil the underlying mechanism is needed.

Discussion

Rumen epithelial microbiota are directly associated with host tissue; thus, it is likely that they play a vital role in mucosal immune systems and ruminant metabolism [12]. Although the last decades have given further insight into the composition of the rumen epithelial microbiota [11], there is limited knowledge of the crosstalk between rumen epithelium-attached bacteria and the host’s cells to date. In this study, we integrated the microbiome, metabolome, and single-cell transcriptome of rumen epithelial tissue to explore the differences in microbiota-host crosstalk between the NB and AD cattle.

In newborn calves, the immune cells were dominated by innate immune cells (dendritic cells), while in adult cattle (as also shown in our multi-tissue study [19]), immune cells were dominated by adaptive immune cells (Th17 cells). This difference between newborns and adults agreed with the results of a previous study that reported Langerhans cells (a dendritic cell type) were present in relative large numbers and T cells occurred at low frequency in the rumen mucosa of the sheep fetus, while T cells (mainly CD4+ T cell subsets) were present in relative large numbers in the adult rumen mucosa [52]. Dendritic cells play an important role in maintaining gastrointestinal homeostasis and priming the adaptive immune system by processing and presenting antigens to T cells [53]. There was a large number of dendritic cells presented in newborn rumen tissues, suggesting that the newborn ruminal environment is primed for T cell recruitment and memory generation. Josefsen et al. [52] inferred that the prevalence of T cells in the adult rumen mucosa may be influenced by antigen leakage through the epithelium. The results of the present study may be attributed to this reason that the rumen of adult cattle had experienced changes in diet (from milk to pellets and hay) and environment, and continuous exposure to feed, environmental antigens, and a variety of symbiotic bacteria may make the rumen the main site of pathogen infection. Previous studies reported that the commensal microbe, segmented filamentous bacterium (SFB), induced the generation of homeostatic Th17 cells that reside at mucosal surfaces where they protect host from pathogenic bacteria in the rodent gastrointestinal tract [54, 55]. However, inconsistent with this, the present study showed a high percentage of Th17 cells in AD rumen compared with NB rumen but no detection of SFB in both groups. This may be due to species specificity of the bacteria-dependent accumulation of Th17 cells. Bifidobacterium adolescentis, the analogously functioning microbes of SFB, that could induce Th17 cells in the murine intestine was found in humans [56]. Therefore, further studies are needed to identify bacterial species from the bovine rumen microbiota capable of inducing Th17 cells.

Single-cell transcriptomics offers the opportunity to compare cell types across species. The rumen epithelium is composed of stratified squamous epithelium, while the human stomach is composed of columnar epithelium [16, 57]. Due to the different anatomical structures of epithelium in the bovine rumen and the human stomach, it is reasonable that bovine rumen epithelial cells are not similar to that of human stomach. A previous study reported the conserved gene expression patterns between skin and rumen based on traditional molecular biology techniques (bulk RNA sequencing) [39]. Although anatomically similar to rumen epithelium, the epithelial cells of human skin (also composed of stratified squamous epithelium) were not similar to those of bovine rumen at single-cell resolution in this study, which suggested that the species difference overrides the cell-type similarity in orthologous gene expression. We observed that the IGFBP3+ EC appeared only in newborn calves considered to be functionally monogastric animals. The IGFBP3 encodes a protein that is known to inhibit the role of IGFs in cell proliferation [58]. This indicates that the disappearance of IGFBP3+ EC may facilitate the rumen epithelium development. Moreover, fibroblasts and vascular endothelial cells were more abundant in the newborn rumen, which may contribute to the rapidly postnatal development of rumen epithelial tissue of newborn calves.

In this study, noticeable changes in the mucosal microbiota were observed between the NB and AD groups, which was reflected by alpha diversity indices (significantly increased in the AD group) and a few genera shared between groups. There were 95 and 141 genera with relative abundance greater than 0.01% in the NB and AD groups, respectively, but only 25 were shared between the two groups, and 10 and 22 genera with relative abundance greater than 1% were found in the NB and AD groups, respectively, but only 1 of them was shared between the two groups. Similar observations have been reported for the goat ruminal epithelial bacteria [15] and the bovine rumen bacterial community [59]. Jiao et al. [15] reported that the alpha diversity indices of goat ruminal epithelial bacteria increased with age after birth and each age group had its distinct epithelial microbiota. Jimi et al. [59] also reported that the diversity of the bacterial community in the bovine rumen fluid increased with age and the similarity of bacteria between 1~3-day-old and other age groups was all very low, along with only a few shared genera. The genera Escherichia-Shigella, Streptococcus, and Bacteroides were dominant in NB rumen tissues, which was similar to the results of a previous study that reported the Escherichia-Shigella, Streptococcus, and Bacteroides were dominant at early after birth (0~28 days), but noticeably decreased at 42 or 70 days after birth in the goat rumen epithelium [15]. The genera Escherichia-Shigella and Streptococcus comprise facultatively anaerobic bacteria, which could create the reduced environment that is required for anaerobic microbes [15, 60, 61]. The predominant genera in adult rumen epithelium mainly included Butyrivibrio, Prevotella, Campylobacter, and Desulfobulbus, which were consistent with the results of previous investigations on rumen mucosa of dairy cows [62]. The Butyrivibrio genus includes the major known cellulolytic species, and the Prevotella genus encompasses a wide array of species that are capable of utilizing different substrates (starches, other non-cellulosic polysaccharides, and simple sugars) [63, 64]. However, little is known about how these genera interacted with the host, which remains to be elucidated.

The most significant transcriptional differences between newborn and adult ruminants were observed in the rumen by comprehensively characterizing transcription of the entire gastrointestinal tract [9]. In this study, compared with the adult rumen, the upregulated DEGs of many epithelial cell subtypes in the newborn rumen were mainly enriched in GO terms related to cell proliferation and development. A previous study found that the higher relative abundance of Megasphaera was involved in promoting rumen epithelial growth of calves [65]. Corresponding to this, the present study also showed a higher relative abundance of Megasphaera in the rumen epithelium of newborn calves compared with AD ones. Megasphaera can convert lactate and glucose into butyrate that is responsible for rumen epithelial proliferation and development [66, 67]. A previous study also suggested that microbes may be involved in the early development of rumen [10], even though the molecular mechanisms behind the development of rumen tissue have been considered ontogenic [68]. These results indicate a potential role of bacteria-driven regulation in the transcriptional activity of rumen epithelial cells; however, further studies are needed to reveal the behind mechanism.

A previous study reported that a co-expressed module gene set involved in “tissue metabolism-related” functions and related to “respiratory electron transport chain” was positively correlated with age of calves and was related to the relative abundance of bacterial genera in the rumen [10]. In the present study, we found that the rumen epithelial cell subtypes had higher metabolism-related functions as well as the oxidation-reduction process at adult stages, which indicates maintaining ROS balance will be very important for the physiological function of cells. A recent study identified a core rumen epithelial microbiota in cattle including the sulfate-reducing bacteria Desulfovibrio [11]. The Desulfovibrio genus converts sulfate to hydrogen sulfide and plays essential roles in reducing ROS [69, 70]. Our study found the Desulfovibrio had a higher relative abundance at the adult stage compared with the newborn stage of the rumen. The relative abundance of Desulfovibrio was positively correlated with concentrations of pyridoxal. Genomic analysis of rumen-associated Desulfovibrio also revealed that it contains genes that are involved in the pyridoxal biosynthesis [71]. The pyridoxal, as an electron donor for sulfate reduction, not only was oxidized by Desulfovibrio to reduce ROS [51] but also prevents host tissue from oxidative stress [72, 73]. Compared with the NB group, adult rumen epithelial cell subtypes showed a higher ability of vitamin B6 (pyridoxal) binding, which suggests that they may uptake more vitamin B6 or participate in B6-mediated signaling pathways. Zhang et al. [74] found that the relative abundance and copy number of Desulfovibrio in the rumen was significantly increased after supplementing antioxidants (e.g., resveratrol) in dairy cattle, suggesting that Desulfovibrio may be involved in the oxidation-reduction process in the rumen. We also noted that the betaine had the highest co-occurrence probability with the Desulfovibrio genus. Previous studies reported that the Desulfovibrio genus can synthesize and uptake betaine to regulate osmotic pressure for adapting to environmental stress (e.g., high salinity) [75, 76]. Betaine, also as a methyl donor, is involved in regulating the remethylation of homocysteine to methionine, and this reaction is catalyzed by betaine-homocysteine methyltransferase (BHMT) [77]. A recent study showed that rumen-protected betaine supplementation can help improve milk yield and milk protein yield in dairy cows, which may be attributed to sparing methionine from being used as a methyl donor and improving the incorporation of methionine into milk protein [78]. However, no significant difference in the expression level of the BHMT gene between cell types of calves and adult cattle was found in this study. Therefore, the mechanisms underlying how the Desulfovibrio interact with the rumen epithelial cell subtypes remain to be elucidated.

In addition to these notable findings, there are some limitations in the current study, which will require additional work and new tools to address. Although previous studies reported that a small sample size should be sufficient for the scRNA-seq [79,80,81], microbiome, and metabolome analysis [15, 82,83,84], other potential associations between the rumen mucosal bacteria taxa and host cells may not be observed due to the limited number of calves (3 animals per group) used in this study. And although we have enabled deeper and more detailed insights into aspects of the rumen epithelium and its microbiota, these associations mainly represent snapshots in time and space. Future studies to collect and analyze samples with larger sample sizes and more time points may be of great importance to completely understand the host-microbe interactions in the rumen.

Conclusions

In summary, we delineated the dynamics of cell-type composition and cellular functions, the relative abundance of the mucosal bacteria, and the potential microbiota-host crosstalk between newborn and adult rumen tissues (Fig. 6). Our integrated analysis of the microbiome, metabolome, and single-cell transcriptome of rumen epithelial tissue provides novel and fundamental insights into the functional completeness of the rumen and may guide future precision feeding in the livestock industry.

Fig. 6
figure 6

Schematic model summarizing the different epithelial microbiota, cell types, and their interactions between NB and AD rumen. NB, newborn; AD, adult

Methods

Animals

All the animal experimental protocols were conducted in compliance with the Zhejiang University Animal Care guidelines and approved by the Animal Care Committee at Zhejiang University (Hangzhou, China). Three adult Holstein dairy cattle with similar age (3.2~4.4 years old), parity (2~3), milk production (18~22 kg/day), and days in milk (181~200 days) were housed in an open barn and fed the same diet. Three newborn calves (1 day after birth) with similar body weight were transferred to a surgery room soon after birth to collect rumen tissue samples. Epithelial microbiome, metabolome, and single-cell transcriptome analysis were conducted in the same animals.

Rumen epithelial tissues acquisition

All animals were humanely euthanized in a surgery room to collect samples from the ventral rumen sac. Samples were rinsed in Dulbecco’s phosphate-buffered saline and quickly transferred into cold magnetic activated cell sorting (MACS) tissue storage solution (Miltenyi Biotec, Bergisch Gladbach, Germany) and sent back to the laboratory immediately for single-cell suspension preparation. Samples used for epithelial microbiome and metabolome analysis were snap-frozen in liquid nitrogen and stored at −80 °C.

Rumen single-cell suspension preparation

After sending back to the laboratory within the cold MACS Tissue Storage Solution, the rumen tissues isolated from adult and newborn cattle were stripped of the outer muscle layers, and then minced into 10 × 0.5mm2 pieces on ice with scissors. The tissues were incubated in a 37°C water bath with 20 mM EDTA for 30 min, and then were rinsed with Dulbecco’s phosphate-buffered saline and chopped into 1 mm pieces. Tissue pieces were transferred to a 15-mL centrifuge tube and re-suspended with 0.25% Trypsin-EDTA (Cat# 25200056, Gibco). After incubated in a 37°C water bath for 5 min, the centrifuge tube containing tissues was inserted into ice for 2 min, and pre-chilled Hank’s balanced salt solution (HBSS) was added to stop the digestion. After centrifuging at 300×g for 2 min at 4°C, the supernatant was discarded, and samples were washed twice with cold HBSS and then re-suspended with dissociation enzymes. Samples were treated with enzymes (1.5mg/ml collagenase I, 1.5ml/ml collagenase IV, 1.5mg/ml dispase, 100U/ml hyaluronidase, and 50U/ml DNase I) for 30 min at 37°C. The digestion was stopped by adding 10% of fetal bovine serum, followed by a filtration step through a 70- and 30-μm SmartStrainer (Miltenyi Biotec, Bergisch Gladbach, Germany). Samples were centrifuged at 300×g for 5 min at 4°C, and then re-suspended in 2 mL of HBSS. Dissociated cells were centrifuged at 300×g for 5 min at 4°C, washed twice with 1×PBS with 0.04% BSA, centrifuged at 300×g for 5 min at 4°C, and re-suspended in 1×PBS with 0.04% BSA.

Viability for single-cell suspension was assessed via trypan blue using a Countess II Automated Cell Counter. If viability for samples were low, the MACS Dead Cell Remove Kit (Miltenyi Biotec, Bergisch Gladbach, Germany) was used to remove dead cells following the manufacturer’s recommendations. Finally, cells were diluted to a concentration of 700–1200 cells/μl with 1× PBS with 0.04% BSA for 10× Genomics sequencing.

cDNA library preparation and sequencing

The 10× Genomics Chromium machine was used for single-cell capture, and the library preparation was performed using the Chromium Single Cell 3′ Reagent Kits v3. A single-cell library was generated separately for each animal. After checking the quality using the Agilent Bioanalyzer High Sensitivity chip, the libraries were sequenced on NovaSeq 6000 platform in a 150-bp pair-ended manner.

Rumen scRNA-seq data analysis

Sequencing results by the NovaSeq 6000 sequencing system were demultiplexed and converted to FASTQ format using Illumina bcl2fastq software. Sample demultiplexing, barcode processing, and single-cell 3′ gene counting was calculated using the Cell Ranger v3.1.0, and scRNA-seq data were aligned to the ARS-UCD1.2 cattle reference genome. The Seurat [85] (version 4.0.3) was used for dimensional reduction, clustering, and analysis of scRNA-seq data. Cells with less than 500 detected genes were considered low quality or empty [22, 86], and cells with more than 4000 detected genes were considered potential doublets [21]. Therefore, for each dataset, we filtered low-quality cells with <500 and >4000 measured genes, UMI counts higher than 50,000, and a mitochondrial gene ratio of higher than 40%. The DoubletFinder [87] package (version 2.0.3) was further used to remove doublets.

For the rumen single-cell atlas, individual newborn and adult Seurat objects were firstly merged into an object. We performed batch correction using Harmony [88] for data integration between samples. The merged object was performed using the NormalizeData, FindVariableFeatures, ScaleData, RunPCA, RunHarmony, FindNeighbors, FindClusters, and RunUMAP functions implemented in the Seurat and Harmony packages. We used the top 30 harmony dimensions and a resolution parameter set to 0.1 for cell clustering. Next, the “FindAllMarkers” function in Seurat was used to determine the marker genes (|‘avg_logFC’| > 0.25 and ‘p_val_adj’ < 0.05) for each cell cluster.

To generate immune cell atlas and non-immune cell atlas in the rumen, we first performed cluster analysis with the top 20 principal components and a resolution of 0.8 for each aforementioned filtered rumen dataset. Cell clusters marked by the immune cell canonical marker gene PTPRC were selected and then used for integration to create an immune cell atlas. In this step, the PTPRC+ cell cluster was not found in the Rumen_NB_1 sample (one of the single-cell datasets of calves), which may be due to the small number of immune cells (only 18 immune cells were found in this sample in the above rumen single-cell atlas). A previous study reported that the sufficiently rare cell type (smaller than about 30 cells) will be algorithmically grouped with a more populous cell type because entire cell clusters, rather than individual cells, were annotated in each tissue [89]. Therefore, there were no cells from the Rumen_NB_1 sample in the immune cell atlas. All individual datasets without PTPRC+ cells were then used for integration to create a non-immune cell atlas. Individual PTPRC+ and PTPRC Seurat objects were merged separately into immune cell and non-immune cell objects, respectively. Each of these merged objects was performed using dimensional reduction, batch-effect correction, and clustering. We performed cluster analysis with the top 50 and 50 harmony dimensions, a resolution of 0.2 and 0.8 for the merged immune and non-immune cell datasets, respectively. Finally, the “FindAllMarkers” function was used to determine the marker genes for each cell cluster. Specifically, at each of the following steps, non-relevant cell types were filtered followed by recalculation of dimensional reduction, batch-effect correction, and clustering. The related marker genes for each cell cluster of immune and non-immune cell atlas are available in the Additional file 2: Table S2 and Table S3, respectively.

Cell-type composition variation analysis

The numbers of cells of each immune cell type in the different groups (NB and AD) were counted and divided by the total number of immune cells in the same group to calculate the percentage of a given immune cell type for each group. The Log2FC between the AD and NB groups was then calculated to identify the cell types altered. The |Log2FC| > 0.5 is considered statistically significant according to the published study [81], which means the cell types were significantly altered between the AD and NB groups.

Cross-species comparison

We downloaded a single-cell RNA-seq dataset of the human stomach from GSE134355 (7 samples) [90] in the Gene Expression Omnibus database and a single-cell RNA-seq dataset of human skin from HRA000395 (3 samples; HRI077736, HRI077737, and HRI077738) [91] in the Genome Sequence Archive database. Firstly, we performed the cell clustering analysis with a resolution parameter setting of 0.8 and 1.0 for the human stomach and human skin, respectively. Next, cell clusters marked by the epithelial cell canonical marker genes were selected in the human stomach and the skin dataset, and the rumen epithelial cell clusters from the non-immune cell atlas were also selected. Finally, the expression matrices (raw counts) of epithelial cell clusters of the two species were merged (cattle rumen and human stomach; cattle rumen and human skin) and orthologous genes were extracted from the data to enable cross-species analysis. The merged expression matrices were used to perform the MetaNeighbour [92] analysis with default parameters. The mean AUROC scores were obtained from MetaNeighbour, and if the AUROC scores between cell types were higher than 0.9, they were considered to be similar or conserved.

Differential gene expression analysis

We used the “FindMarkers” function in Seurat to identify DEGs for each cell type between NB and AD groups. The normalized read count that was obtained by executing “NormalizeData” and “ScaleData” functions in Seurat was used for the differential gene expression analysis. And the log fold change (Log2FC) and p-value of each DEG were calculated by using the Wilcoxon rank-sum test as implemented in the “FindMarkers” function. The p-value adjustment was performed using Bonferroni correction based on the total number of genes in the dataset. The cell types with fewer than 50 cells in the NB or AD groups were filtered before executing the differential gene expression analysis. Only genes expressed in more than 15% of the cells in the specific cell type were considered. DEGs between the NB and AD groups were identified to generate upregulated and downregulated DEG datasets (|logFC| > 0.5, adjusted p-value < 0.05) for each cell type.

GO term analysis

GO term enrichment analysis was performed using the function enrichGO in clusterProfiler R package [93] based on the dataset “org.Bt.eg.db.” Dot plot of representative GO terms based on the upregulated and downregulated DEGs of cell type were generated with ggplot2 [94].

16S rRNA gene sequencing and analysis

Total DNA of the epithelial microbial community was extracted from each rumen epithelial tissue using the E.Z.N.A.® DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s instructions. The qualities and quantities of the DNA samples were determined with a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA). The hypervariable region V3-V4 of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCRs were performed with the following program: 95°C for 3 min; 28 cycles (for region V3-V4 of 16S rRNA gene); followed by 72°C for 10 min. PCR products were extracted from 2% agarose gel and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instruction and quantified using Quantus™ Fluorometer (Promega, USA). Amplicon sequencing was performed on an Illumina MiSeq platform (Illumina, San Diego, USA) using the paired-end 2×300-bp protocol.

The paired-end reads were merged using FLASH (version 1.2.11) [95]. The 16S rRNA gene sequencing analysis was performed using the QIIME2 [96]. Reads were truncated at the first instance of a quality score less than 20, and then reads with a length less than 50bp or containing N base were removed. In the DADA2 plugin [97], the “filterAndTrim” function with the parameters “MaxEE” setting to 2 and the “truncQ” setting to 0 was used to perform further quality control and chimera was removed using the “consensus” method of the “removeBimeraDenovo” function with the default parameters, and finally, an ASV feature table was produced. The ASV feature table was used for taxonomic identification using a Naiv̈e Bayes classifier trained on the Silva database (Release138, http://www.arb-silva.de) and clustering at 99%. The Shannon and Chao1 indexes of ASVs were obtained in QIIME2. Linear discriminant analysis (LDA) effect size (LEfSe) was used to identify bacterial taxa that were significantly (LDA > 3.5, p < 0.05) enriched in NB or AD rumen tissues. The data used for the bacterial composition analysis and the LEfSe analysis were based on the ASVs that were merged to the genus level.

Metabolomics analysis of rumen epithelial tissues

In the current study, the metabolomics analysis experiment on the newborn rumen tissue was performed at the same time as and using the same batch of reagents as that for the adult rumen tissue in our other study [19]. Specifically, rumen epithelial tissues were homogenized with 1000 μl of ice-cold methanol/water (70%, v/v) and cold steel balls for 3 min at 30 Hz. The tissues were whirled for 1 min without steel balls, and then centrifuged at 4°C, 12,000 rpm for 10 min after 15 min standing. The supernatant was collected for LC-MS/MS analysis.

To identify and quantify as many metabolites as possible in rumen tissues, firstly, all sample extracts are mixed in equal amounts for LC-QTOF-MS/MS experiment to detect metabolites based on the MWDB database built by standard material and public database MHK (including Metlin, HMDB, KEGG database information, secondary spectrum, retention time). Then the information of identified metabolites including multiple ion pair information and retention time is combined with the MWDB database and the accurate quantification information of samples was obtained by Q-Trap. Next, ESI-Q TRAP-MS/MS analysis was performed using the supernatant of each sample. The ESI source operation parameters were as follows: source temperature 500°C; ion spray voltage 5500 V (positive), −4500 V (negative); ion source gas I, gas II, curtain gas were set at 50, 50, and 25 psi, respectively; the collision gas was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/l polypropylene glycol solutions in QQQ and LIT modes, respectively. A specific set of MRM transitions were monitored for each period according to the metabolites eluted within this period.

The OPLS-DA (orthogonal projections to latent structures discriminant analysis) was performed for the identification of the significantly different metabolites between NB and AD groups. Different metabolites were determined by VIP ≥ 1 and absolute of Log2FC ≥ 1. In order to avoid overfitting, a permutation test (200 permutations) was performed.

Microbe–metabolite vectors (mmvec) neural network analysis

To predict the co-occurrence probabilities between microbes at the genus level and metabolites, the mmvec neural networks analysis was applied [47]. According to the analysis pipeline, the “paired-omics” function with the parameters “epochs” setting to 100 and the “learning-rate” setting to 1e−3 was performed to estimate the conditional probability that each metabolite is present given the presence of specific microbes based on the microbial sequence counts and the metabolite relative concentration. Only microbes at the genus level that appear in at least 2 samples were included in the mmvec neural networks analysis.

Gene set scoring analysis

The genes of the “vitamin B6 binding” gene set are listed in Additional file 2: Table S9. The signature score of each gene set in each cell type from the NB or AD group was calculated using the AddModuleScore function in the Seurat R package. The differences in the signature scores between two groups of cells were evaluated by a two-sided Wilcoxon rank-sum test.

Fluorescence in situ hybridization

The rumen tissues collected from the ventral rumen sac of newborn and adult cattle were fixed in a 4% paraformaldehyde fix solution (Beyotime, Shanghai, China) for 24 h and then embedded in paraffin for sectioning. The slices were deparaffinized, dehydrated, and then boiled in the antigen retrieval solution (Cat#G1202, Servicebio) for 10 min. After cooling, slices were added with proteinase K solution (Cat#G1234, Servicebio) and incubated at 37°C for 20 min. Slices were washed three times with phosphate-buffered saline and then added with hybridization buffer (20% formamide) and incubated at 37°C for 1 h. For bacterial FISH, after removing the hybridization buffer, slices were incubated with hybridization buffer (20% formamide) containing the bacterial probe (1μM) overnight at 37°C, and then slices were washed using 2×Saline sodium citrate (SSC) for 10 min, 1×SSC twice (5 min each time) at 37°C, and 0.5×SSC for 10 min at room temperature, respectively. Nuclei were counterstained with DAPI for 8min in the dark. Images were taken using a Pannoramic DESK scanning microscope and analyzed by Caseviewer (version 2.3). The bacterial probe specific to Desulfovibrio used for this study was 5′-FAM-GGTCGCCCCCCGACACCT -FAM-3′ (Cat# SPZL000913, EXONBIO).

Quantification of bacterial density

Absolute quantitative real-time PCR was performed to determine the bacterial densities for the genera Enterococcus, Butyrivibrio, Prevotella, and Desulfovibrio between newborn and adult rumen epithelium. Primers are specific to the genera Enterococcus: forward primer (CCCTTATTGTTAGTTGCCATCATT) and reverse primer (ACTCGTTGTACTTCCCATTGT) [98]. Primers specific to the genera Butyrivibrio: forward primer (GYGAAGAAGTATTTCGGTAT) and reverse primer (CCAACACCTAGTATTCATC) [99]. Primers are specific to the genera Prevotella: forward primer (GGTTCTGAGAGGAAGGTCCCC) and reverse primer (TCCTGCACGCTACTTGGCTG) [12]. Primers specific to the genera Desulfovibrio: forward primer (ACCTGCTGGAACTGCAARA (R = G or A)) and reverse primer (GTGGAAGCCCACGCTGTT) [100].

Statistical analysis

The experimental data were performed using a two-tailed Student’s test in GraphPad Prism software (version 8.0) to compare the differences between NB and AD groups assuming equal variance. P < 0.05 is considered statistically significant.

Availability of data and materials

The single-cell RNA sequencing raw data of the adult rumen were previously analyzed in our cross-tissue single-cell transcriptomic landscape study [19] and are available under accession number SRP321626 [101] in the NCBI Sequence Read Archive (SRA) database. The raw single-cell RNA sequencing data of the newborn cattle rumen, and the merged processing files of single-cell RNA sequencing data from the newborn and adult cattle rumen tissues, have been deposited to the Gene Expression Omnibus database under accession number GSE183285 [102]. The raw metabolomic data from the newborn and adult rumen have been deposited to the Open Archive for Miscellaneous Data of National Genomics Data Center (NGDC) under accession number OMIX002037 [103] and also been deposited in the MetaboLights database with the identifier MTBLS6501 [104]. The raw files of the 16S rRNA gene sequencing data of the newborn and adult rumen have been deposited to the SRA database under accession number PRJNA846365 [105]. The processed sequencing data on the expression profiles (CPM) of the fibroblast marker genes in newborn and adult cattle rumen were collected from the GSE74329 [106] and GSE78197 [107], respectively. Other relevant data are available upon request. This study did not generate any unique code.

Abbreviations

AD:

Adult

BC:

Basal cells

BC1:

Basal cell 1

BC2:

Basal cell 2

BC3:

Basal cell 3

cg-like SC:

Channel-gap-like spinous cells

CPM:

Counts per million

DC:

Dendritic cells

DEGs:

Differentially expressed genes

EC:

Epithelial cells

FISH:

Fluorescence in situ hybridization

GC:

Granule cells

GC1:

Granule cell 1

GC2:

Granule cell 2

GC3:

Granule cell 3

GO:

Gene Ontology

HBSS:

Hank’s balanced salt solution

MACS:

Magnetic activated cell sorting

MC:

Mitotic cell

MC1:

Mitotic cell 1

MC2:

Mitotic cell 2

mmvec:

Microbe–metabolite vectors

NB:

Newborn

NKT:

Natural killer T cells

ROS:

Reactive oxygen species

SC:

Spinous cells

SC1:

Spinous cell 1

SC2:

Spinous cell 2

SFB:

Segmented filamentous bacterium

scRNA-seq:

Single-cell RNA sequencing

VEC:

Vascular endothelial cells

VIP:

Variable importance in projection.

References

  1. Vigne JD. The origins of animal domestication and husbandry: a major change in the history of humanity and the biosphere. C R Biol. 2011;334(3):171–81.

    Article  Google Scholar 

  2. Connor EE, Li RW, Baldwin RL, Li C. Gene expression in the digestive tissues of ruminants and their relationships with feeding and digestive processes. Animal. 2010;4(7):993–1007.

    Article  CAS  Google Scholar 

  3. Matthews C, Crispie F, Lewis E, Reid M, O'Toole PW, Cotter PD. The rumen microbiome: a crucial consideration when optimising milk and meat production and nitrogen utilisation efficiency. Gut Microbes. 2019;10(2):115–32.

    Article  CAS  Google Scholar 

  4. Hofmann RR. Evolutionary steps of ecophysiological adaptation and diversification of ruminants: a comparative view of their digestive system. Oecologia. 1989;78(4):443–57.

    Article  CAS  Google Scholar 

  5. Svartström O, Alneberg J, Terrapon N, Lombard V, de Bruijn I, Malmsten J, et al. Ninety-nine de novo assembled genomes from the moose (Alces alces) rumen microbiome provide new insights into microbial plant biomass degradation. ISME J. 2017;11(11):2538–51.

    Article  Google Scholar 

  6. Heinrichs A, Lesmeister K. Rumen development in the dairy calf. Calf and heifer rearing: principles of rearing the modern dairy heifer from calf to calving 60th. Nottingham, UK: University of Nottingham Easter School in Agricultural Science; 2005. p. 53–65.

    Google Scholar 

  7. Braun U, Kruger S, Hassig M. Ultrasonographic examination of the reticulum, rumen, omasum and abomasum during the first 100 days of life in calves. Res Vet Sci. 2013;95(2):326–33.

    Article  CAS  Google Scholar 

  8. Stobo I, Roy J, Gaston HJ. Rumen development in the calf: 1. The effect of diets containing different proportions of concentrates to hay on rumen development. Br J Nutr. 1966;20(2):171–88.

    Article  CAS  Google Scholar 

  9. Bush SJ, Mc Culloch MEB, Muriuki C, Salavati M, Davis GM, Farquhar IL, et al. Comprehensive Transcriptional Profiling of the Gastrointestinal Tract of Ruminants from Birth to Adulthood Reveals Strong Developmental Stage Specific Gene Expression. G3. 2019;9(2):359–73.

    Article  CAS  Google Scholar 

  10. Malmuthuge N, Liang G, Guan LL. Regulation of rumen development in neonatal ruminants through microbial metagenomes and host transcriptomes. Genome Biol. 2019;20(1):172.

    Article  Google Scholar 

  11. Anderson CJ, Koester LR, Schmitz-Esser S. Rumen epithelial communities share a core bacterial microbiota: a meta-analysis of 16S rRNA gene Illumina MiSeq sequencing datasets. Front Microbiol. 2021;12:625400.

    Article  Google Scholar 

  12. Malmuthuge N, Griebel PJ, Guan LL. Taxonomic identification of commensal bacteria associated with the mucosa and digesta throughout the gastrointestinal tracts of preweaned calves. Appl Environ Microbiol. 2014;80(6):2021–8.

    Article  Google Scholar 

  13. McCowan RP, Cheng KJ, Bailey CB, Costerton JW. Adhesion of bacteria to epithelial cell surfaces within the reticulo-rumen of cattle. Appl Environ Microbiol. 1978;35(1):149–55.

    Article  CAS  Google Scholar 

  14. Rieu F, Fonty G, Gaillard B, Gouet P. Electron microscopy study of the bacteria adherent to the rumen wall in young conventional lambs. Can J Microbiol. 1990;36(2):140–4.

    Article  CAS  Google Scholar 

  15. Jiao J, Huang J, Zhou C, Tan Z. Taxonomic identification of ruminal epithelial bacterial diversity during rumen development in goats. Appl Environ Microbiol. 2015;81(10):3502–9.

    Article  CAS  Google Scholar 

  16. Steele MA, Penner GB, Chaucheyras-Durand F. Development and physiology of the rumen and the lower gut: Targets for improving gut health. J Dairy Sci. 2016;99(6):4955–66.

    Article  CAS  Google Scholar 

  17. Haber AL, Biton M, Rogel N, Herbst RH, Shekhar K, Smillie C, et al. A single-cell survey of the small intestinal epithelium. Nature. 2017;551(7680):333–9.

    Article  CAS  Google Scholar 

  18. Gao Y, Fang L, Baldwin R.L. VI, Connor EE, Cole JB, Tassell CPV, et al. Single-cell transcriptomic analyses of dairy cattle ruminal epithelial cells during weaning. Genomics. 2021;113(4):2045–55.

  19. Wu JJ, Zhu S, Gu F, Valencak TG, Liu JX, Sun HZ. Cross-tissue single-cell transcriptomic landscape reveals the key cell subtypes and their potential roles in the nutrient absorption and metabolism in dairy cattle. J Adv Res. 2022;37:1–18.

    Article  CAS  Google Scholar 

  20. Xue M-Y, Wu J-J, Xie Y-Y, Zhu S-L, Zhong Y-F, Liu J-X, et al. Investigation of fiber utilization in the rumen of dairy cows based on metagenome-assembled genomes and single-cell RNA sequencing. Microbiome. 2022;10(1):11.

    Article  CAS  Google Scholar 

  21. Buechler MB, Pradhan RN, Krishnamurty AT, Cox C, Calviello AK, Wang AW, et al. Cross-tissue organization of the fibroblast lineage. Nature. 2021;593(7860):575–9.

    Article  CAS  Google Scholar 

  22. Kalucka J, de Rooij L, Goveia J, Rohlenova K, Dumas SJ, Meta E, et al. Single-cell transcriptome atlas of murine endothelial cells. Cell. 2020;180(4):764–779 e720.

    Article  CAS  Google Scholar 

  23. Han X, Zhou Z, Fei L, Sun H, Wang R, Chen Y, et al. Construction of a human cell landscape at single-cell level. Nature. 2020;581(7808):303–9.

    Article  CAS  Google Scholar 

  24. Cunha P, Le Vern Y, Gitton C, Germon P, Foucras G, Rainard P. Expansion, isolation and first characterization of bovine Th17 lymphocytes. Sci Rep. 2019;9(1):1–14.

    Article  Google Scholar 

  25. Pizzolato G, Kaminski H, Tosolini M, Franchini DM, Pont F, Martins F, et al. Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVdelta1 and TCRVdelta2 gammadelta T lymphocytes. Proc Natl Acad Sci U S A. 2019;116(24):11906–15.

    Article  CAS  Google Scholar 

  26. Han X, Wang R, Zhou Y, Fei L, Sun H, Lai S, et al. Mapping the Mouse Cell Atlas by Microwell-Seq. Cell. 2018;172(5):1091–1107 e1017.

    Article  CAS  Google Scholar 

  27. Young MD, Mitchell TJ, Vieira Braga FA, Tran MGB, Stewart BJ, Ferdinand JR, et al. Single-cell transcriptomes from human kidneys reveal the cellular identity of renal tumors. Science. 2018;361(6402):594–9.

    Article  CAS  Google Scholar 

  28. Talker SC, Baumann A, Barut GT, Keller I, Bruggmann R, Summerfield A. Precise delineation and transcriptional characterization of bovine blood dendritic-cell and monocyte subsets. Front Immunol. 2018;9:2505.

    Article  Google Scholar 

  29. Villani A-C, Satija R, Reynolds G, Sarkizova S, Shekhar K, Fletcher J, et al. Single-cell RNA-seq reveals new types of human blood dendritic cells, monocytes, and progenitors. Science. 2017;356(6335).

  30. Aschenbach JR, Zebeli Q, Patra AK, Greco G, Amasheh S, Penner GB. Symposium review: The importance of the ruminal epithelial barrier for a healthy and productive cow. J Dairy Sci. 2019;102(2):1866–82.

    Article  CAS  Google Scholar 

  31. Graham C, Simmons NL. Functional organization of the bovine rumen epithelium. Am J Physiol Regul Integr Comp Physiol. 2005;288(1):R173–81.

    Article  CAS  Google Scholar 

  32. Haensel D, Jin S, Sun P, Cinco R, Dragan M, Nguyen Q, et al. Defining epidermal basal cell states during skin homeostasis and wound healing using single-cell transcriptomics. Cell Rep. 2020;30(11):3932–47 e3936.

    Article  CAS  Google Scholar 

  33. Villemin J-P, Lorenzi C, Cabrillac M-S, Oldfield A, Ritchie W, Luco RF. A cell-to-patient machine learning transfer approach uncovers novel basal-like breast cancer prognostic markers amongst alternative splice variants. BMC Biol. 2021;19(1):70.

    Article  CAS  Google Scholar 

  34. Yohe TT, Tucker HLM, Parsons CLM, Geiger AJ, Akers RM, Daniels KM. Short communication: initial evidence supporting existence of potential rumen epidermal stem and progenitor cells. J Dairy Sci. 2016;99(9):7654–60.

    Article  CAS  Google Scholar 

  35. He H, Suryawanshi H, Morozov P, Gay-Mimbrera J, Del Duca E, Kim HJ, et al. Single-cell transcriptome analysis of human skin identifies novel fibroblast subpopulation and enrichment of immune subsets in atopic dermatitis. J Allergy Clin Immunol. 2020;145(6):1615–28.

    Article  CAS  Google Scholar 

  36. Robitaille H, Proulx R, Robitaille K, Blouin R, Germain L. The mitogen-activated protein kinase kinase kinase dual leucine zipper-bearing kinase (DLK) acts as a key regulator of keratinocyte terminal differentiation. J Biol Chem. 2005;280(13):12732–41.

    Article  CAS  Google Scholar 

  37. Zou Z, Long X, Zhao Q, Zheng Y, Song M, Ma S, et al. A single-cell transcriptomic atlas of human skin aging. Dev Cell. 2020;56(3):383–397.e8.

    Article  Google Scholar 

  38. Goodenough DA, Paul DL. Gap junctions. Cold Spring Harb Perspect Biol. 2009;1(1):a002576.

    Article  Google Scholar 

  39. Xiang R, McNally J, Rowe S, Jonker A, Pinares-Patino CS, Oddy VH, et al. Gene network analysis identifies rumen epithelial cell proliferation, differentiation and metabolic pathways perturbed by diet and correlated with methane production. Sci Rep. 2016;6:39022.

    Article  CAS  Google Scholar 

  40. Wang B, Wang D, Wu X, Cai J, Liu M, Huang X, et al. Effects of dietary physical or nutritional factors on morphology of rumen papillae and transcriptome changes in lactating dairy cows based on three different forage-based diets. BMC Genomics. 2017;18(1):353.

    Article  Google Scholar 

  41. Liu J, Zhang ML, Zhang RY, Zhu WY, Mao SY. Comparative studies of the composition of bacterial microbiota associated with the ruminal content, ruminal epithelium and in the faeces of lactating dairy cows. J Microbial Biotechnol. 2016;9(2):257–68.

    Article  CAS  Google Scholar 

  42. Guo W, van Niekerk J, Zhou M, Steele M, Guan LL. Longitudinal assessment revealed the shifts in rumen and colon mucosal-attached microbiota of dairy calves during weaning transition. J Dairy Sci. 2021;104(5):5948–63.

    Article  CAS  Google Scholar 

  43. Seshadri R, Leahy SC, Attwood GT, Teh KH, Lambie SC, Cookson AL, et al. Cultivation and sequencing of rumen microbiome members from the Hungate1000 Collection. Nat Biotechnol. 2018;36(4):359–67.

    Article  CAS  Google Scholar 

  44. Newbold CJ, Ramos-Morales E. Review: ruminal microbiome and microbial metabolome: effects of diet and ruminant host. Animal. 2020;14:s78–86.

    Article  CAS  Google Scholar 

  45. Liu K, Zhang Y, Yu Z, Xu Q, Zheng N, Zhao S, et al. Ruminal microbiota–host interaction and its effect on nutrient metabolism. Anim Nutr. 2021;7(1):49–55.

    Article  CAS  Google Scholar 

  46. Han S, Van Treuren W, Fischer CR, Merrill BD, DeFelice BC, Sanchez JM, et al. A metabolomics pipeline for the mechanistic interrogation of the gut microbiome. Nature. 2021;595(7867):415–20.

    Article  CAS  Google Scholar 

  47. Morton JT, Aksenov AA, Nothias LF, Foulds JR, Quinn RA, Badri MH, et al. Learning representations of microbe-metabolite interactions. Nat Methods. 2019;16(12):1306–14.

    Article  CAS  Google Scholar 

  48. Eliot AC, Kirsch JF. Pyridoxal phosphate enzymes: mechanistic, structural, and evolutionary considerations. Annu Rev Biochem. 2004;73(1):383–415.

    Article  CAS  Google Scholar 

  49. Li J, Yin L, Wang L, Li J, Huang P, Yang H, et al. Effects of vitamin B6 on growth, diarrhea rate, intestinal morphology, function, and inflammatory factors expression in a high-protein diet fed to weaned piglets. J Anim Sci. 2019;97(12):4865–74.

    Article  Google Scholar 

  50. Komatsu S-i, Watanabe H, Oka T, Tsuge H, Nii H, Kato N. Vitamin B-6–supplemented diets compared with a low vitamin B-6 diet suppress azoxymethane-induced colon tumorigenesis in mice by reducing cell proliferation. J Nutr. 2001;131(8):2204–7.

    Article  CAS  Google Scholar 

  51. Zellner G, Kneifel H, Winter J. Oxidation of benzaldehydes to benzoic acid derivatives by three Desulfovibrio strains. Appl Environ Microbiol. 1990;56(7):2228–33.

    Article  CAS  Google Scholar 

  52. Josefsen TD, Landsverk T. T cell subsets and Langerhans cells in the forestomach mucosa of adult sheep and sheep foetuses. Vet Immunol Immunopathol. 1996;51(1-2):101–11.

    Article  CAS  Google Scholar 

  53. Ruane DT, Lavelle EC. The role of CD103+ dendritic cells in the intestinal mucosal immune system. Front Immunol. 2011;2:25.

    Article  Google Scholar 

  54. Ivanov II, Atarashi K, Manel N, Brodie EL, Shima T, Karaoz U, et al. Induction of intestinal Th17 cells by segmented filamentous bacteria. Cell. 2009;139(3):485–98.

    Article  CAS  Google Scholar 

  55. Kumar P, Monin L, Castillo P, Elsegeiny W, Horne W, Eddens T, et al. Intestinal interleukin-17 receptor signaling mediates reciprocal control of the gut microbiota and autoimmune inflammation. Immunity. 2016;44(3):659–71.

    Article  CAS  Google Scholar 

  56. Tan TG, Sefik E, Geva-Zatorsky N, Kua L, Naskar D, Teng F, et al. Identifying species of symbiont bacteria from the human gut that, alone, can induce intestinal Th17 cells in mice. Proc Natl Acad Sci U S A. 2016;113(50):E8141–50.

    Article  CAS  Google Scholar 

  57. Han S, Fink J, Jörg DJ, Lee E, Yum MK, Chatzeli L, et al. Defining the identity and dynamics of adult gastric isthmus stem cells. Cell Stem Cell. 2019;25(3):342–356.e347.

    Article  CAS  Google Scholar 

  58. Nishihara K, Suzuki Y, Kim D, Roh S. Growth of rumen papillae in weaned calves is associated with lower expression of insulin-like growth factor-binding proteins 2, 3, and 6. Anim Sci J. 2019;90(9):1287–92.

    Article  CAS  Google Scholar 

  59. Jami E, Israel A, Kotser A, Mizrahi I. Exploring the bovine rumen bacterial community from birth to adulthood. ISME J. 2013;7(6):1069–79.

    Article  Google Scholar 

  60. Chu J, Zhang T, He K. Cariogenicity features of Streptococcus mutans in presence of rubusoside. BMC Oral Health. 2016;16(1):54.

    Article  Google Scholar 

  61. Jost T, Lacroix C, Braegger CP, Chassard C. New insights in gut microbiota establishment in healthy breast fed neonates. PLoS One. 2012;7(8):e44595.

    Article  CAS  Google Scholar 

  62. Mao S, Zhang M, Liu J, Zhu W. Characterising the bacterial microbiota across the gastrointestinal tracts of dairy cattle: membership and potential function. Sci Rep. 2015;5:16116.

    Article  CAS  Google Scholar 

  63. Mizrahi I, Jami E. Review: The compositional variation of the rumen microbiome and its effect on host performance and methane emission. Animal. 2018;12(s2):s220–32.

    Article  CAS  Google Scholar 

  64. Stevenson DM, Weimer PJ. Dominance of Prevotella and low abundance of classical ruminal bacterial species in the bovine rumen revealed by relative quantification real-time PCR. Appl Microbiol Biotechnol. 2007;75(1):165–74.

    Article  CAS  Google Scholar 

  65. Lin L, Xie F, Sun D, Liu J, Zhu W, Mao S. Ruminal microbiome-host crosstalk stimulates the development of the ruminal epithelium in a lamb model. Microbiome. 2019;7(1):83.

    Article  Google Scholar 

  66. Diao Q, Zhang R, Fu T. Review of strategies to promote rumen development in calves. Animals. 2019;9(8).

  67. Muya M, Nherera F, Miller K, Aperce C, Moshidi P, Erasmus LJ. Effect of Megasphaera elsdenii NCIMB 41125 dosing on rumen development, volatile fatty acid production and blood β-hydroxybutyrate in neonatal dairy calves. J Anim Physiol Anim Nutr. 2015;99(5):913–8.

    Article  CAS  Google Scholar 

  68. Vi RB, McLeod K, Klotz J, Heitmann R. Rumen development, intestinal growth and hepatic metabolism in the pre-and postweaning ruminant. J Dairy Sci. 2004;87:E55–65.

    Article  Google Scholar 

  69. Amann R, Stromley J, Devereux R, Key R, Stahl D. Molecular and microscopic identification of sulfate-reducing bacteria in multispecies biofilms. Appl Environ Microbiol. 1992;58(2):614–23.

    Article  CAS  Google Scholar 

  70. Wallace JL. Physiological and pathophysiological roles of hydrogen sulfide in the gastrointestinal tract. Antioxid Redox Signal. 2010;12(9):1125–33.

    Article  CAS  Google Scholar 

  71. Campbell AG, Campbell JH, Schwientek P, Woyke T, Sczyrba A, Allman S, et al. Multiple single-cell genomes provide insight into functions of uncultured deltaproteobacteria in the human oral cavity. PLoS One. 2013;8(3):e59361.

    Article  CAS  Google Scholar 

  72. Shen J, Lai C-Q, Mattei J, Ordovas JM, Tucker KL. Association of vitamin B-6 status with inflammation, oxidative stress, and chronic inflammatory conditions: the Boston Puerto Rican Health Study. Am J Clin Nutr. 2009;91(2):337–42.

    Article  Google Scholar 

  73. Anand SS. Protective effect of vitamin B6 in chromium-induced oxidative stress in liver. J Appl Toxicol. 2005;25(5):440–3.

    Article  CAS  Google Scholar 

  74. Zhang R, Zhang W, Bi Y, Tu Y, Ma T, Dong L, et al. Sanguinarine and resveratrol affected rumen fermentation parameters and bacterial community in calves. Anim Feed Sci Technol. 2019;251:64–75.

    Article  CAS  Google Scholar 

  75. Welsh DT, Lindsay YE, Caumette P, Herbert RA, Hannan J. Identification of trehalose and glycine betaine as compatible solutes in the moderately halophilic sulfate reducing bacterium, Desulfovibrio halophilus. FEMS Microbiol Lett. 1996;140(2-3).

  76. Ollivier B, Caumette P, Garcia JL, Mah RA. Anaerobic bacteria from hypersaline environments. Microbiol Rev. 1994;58(1):27–38.

    Article  CAS  Google Scholar 

  77. Kidd MT, Ferket PR, Garlich JD. Nutritional and osmoregulatory functions of betaine. World Poultry Sci J. 1997;53(2):125–39.

    Article  Google Scholar 

  78. Wang C, Liu H, Wang C, Liu J, Liu H. Effects of dietary rumen-protected betaine on lactation performance and serum metabolites of mid-lactation Holstein dairy cows. J Agric Food Chem. 2020;68(46):13154–9.

    Article  CAS  Google Scholar 

  79. Guo J, Grow EJ, Mlcochova H, Maher GJ, Lindskog C, Nie X, et al. The adult human testis transcriptional cell atlas. Cell Res. 2018;28(12):1141–57.

    Article  CAS  Google Scholar 

  80. Wang Y, Song W, Wang J, Wang T, Xiong X, Qi Z, et al. Single-cell transcriptome analysis reveals differential nutrient absorption functions in human intestine. J Exp Med. 2020;217(2):e20191130.

    Article  Google Scholar 

  81. Ma S, Sun S, Geng L, Song M, Wang W, Ye Y, et al. Caloric restriction reprograms the single-cell transcriptional landscape of Rattus norvegicus aging. Cell. 2020;180(5):984–1001 e1022.

    Article  CAS  Google Scholar 

  82. Yeoman CJ, Ishaq SL, Bichi E, Olivo SK, Lowe J, Aldridge BM. Biogeographical differences in the influence of maternal microbial sources on the early successional development of the bovine neonatal gastrointestinal tract. Sci Rep. 2018;8(1):3197.

    Article  Google Scholar 

  83. Zhao Y, Shang Y, Ren Y, Bie Y, Qiu Y, Yuan Y, et al. Omics study reveals abnormal alterations of breastmilk proteins and metabolites in puerperant women with COVID-19. Signal Transduct Target Ther. 2020;5(1):247.

    Article  CAS  Google Scholar 

  84. Ding Z, Wang W, Zhang K, Ming F, Yangdai T, Xu T, et al. Novel scheme for non-invasive gut bioinformation acquisition with a magnetically controlled sampling capsule endoscope. Gut. 2021;70(12):2297–306.

    Article  CAS  Google Scholar 

  85. Hao Y, Hao S, Andersen-Nissen E, Mauck WM III, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573–3587.e29.

    Article  CAS  Google Scholar 

  86. Travaglini KJ, Nabhan AN, Penland L, Sinha R, Gillich A, Sit RV, et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature. 2020;587(7835):619–25.

    Article  CAS  Google Scholar 

  87. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8(4):329–337 e324.

    Article  CAS  Google Scholar 

  88. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–96.

    Article  CAS  Google Scholar 

  89. Tabula Muris C, Overall C, Logistical C, Organ C, processing, Library p, sequencing, Computational data a, Cell type a, Writing g et al., Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature. 2018; 562(7727):367-372.

  90. Guo GHX, Zhou Z, Fei L, Sun H, Wang R, Wang J, et al. Construction of a human cell landscape by single-cell mRNA-seq. Gene Expression Omnibus https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134355. 2020.

  91. Liu GH. A single-cell transcriptomic atlas of human skin aging. Genome Sequence Archive https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000395. 2020.

  92. Crow M, Paul A, Ballouz S, Huang ZJ, Gillis J. Characterizing the replicability of cell types defined by single cell RNA-sequencing data using MetaNeighbor. Nat Commun. 2018;9(1):884.

    Article  Google Scholar 

  93. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.

    Article  CAS  Google Scholar 

  94. Wickham H. ggplot2: elegant graphics for data analysis: springer; 2016.

    Book  Google Scholar 

  95. Magoč T, Salzberg SL. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27(21):2957–63.

    Article  Google Scholar 

  96. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37(8):852–7.

    Article  CAS  Google Scholar 

  97. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13(7):581–3.

    Article  CAS  Google Scholar 

  98. Liu X, Cheng Y, Shao L, Ling Z. Alterations of the predominant fecal microbiota and disruption of the gut mucosal barrier in patients with early-stage colorectal cancer. Biomed Res Int. 2020;2020:2948282.

    Google Scholar 

  99. Boeckaert C, Vlaeminck B, Fievez V, Maignien L, Dijkstra J, Boon N. Accumulation of trans C18:1 fatty acids in the rumen after dietary algal supplementation is associated with changes in the Butyrivibrio community. Appl Environ Microbiol. 2008;74(22):6923–30.

    Article  CAS  Google Scholar 

  100. Giacomucci L, Purdy KJ, Zanardini E, Polo A, Cappitelli F. A new non-degenerate primer pair for the specific detection of the nitrite reductase gene nrfA in the genus Desulfovibrio. J Mol Microbiol Biotechnol. 2012;22(6):345–51.

    CAS  Google Scholar 

  101. Wu J, Zhu S, Sun H. Raw single-cell RNA sequencing data of adult rumen tissues. Sequence Read Archive https://www.ncbi.nlm.nih.gov/sra/?term=SRP321626. 2021.

  102. Wu J, Zhu S, Gu F, Sun H. Microbiota-host crosstalk mechanism in the newborn and adult rumen at single-cell resolution. Gene Expression Omnibus, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE183285. 2022.

  103. Gu F. Microbiota-host crosstalk mechanism in the newborn and adult rumen at single-cell resolution. Open Archive for Miscellaneous Data of National Genomics Data Center (NGDC), https://ngdc.cncb.ac.cn/omix/release/OMIX002037. 2022.

  104. Wu J. MTBLS6501: Microbiota-host crosstalk mechanism in the newborn and adult rumen at single-cell resolution. Metabolights. https://www.ebi.ac.uk/metabolights/MTBLS6501. 2022.

  105. Zhejiang University. Bovine gut 16s amplicon Raw sequence reads. Sequence Read Archive, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA846365/. 2022.

  106. Liang G. Transcriptome anlaysis of gastrointestinal tract of pre-weaned cavles. Gene Expression Omnibus https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74329. 2016.

  107. Guan LL, Wang D. RNA sequencing reveals differential expression of genes associated with an altered morphology of rumen papillae in lactating dairy cows fed diets with various forage sources. Gene Expression Omnibus https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78197. 2017.

Download references

Acknowledgements

We thank LC Sciences for deep sequencing experiments. The authors would also like to thank the Core Facilities, Zhejiang University School of Medicine for technical support. The authors would also like to thank Ming-Yuan Xue, Yi-Fan Zhong, and Lu-Yi Jiang of the Institute of Dairy Science, College of Animal Sciences of Zhejiang University for their assistance in sampling.

Funding

This project was supported by the National Natural Science Foundation of China (32002207, 31872380, Beijing), the China Agriculture (Dairy) Research System (CARS-36, Beijing), the Fundamental Research Funds for the Zhejiang Provincial Universities (2021XZZX027), and the “Hundred Talents Program” Research Professor Start-up Fund of Zhejiang University.

Author information

Authors and Affiliations

Authors

Contributions

JJW and HZS designed the study. JJW, HZS, SLZ, and FFG performed sampling. JJW, SLZ, and HZS analyzed the data. JJW and YFT conducted validation experiments. JJW wrote and revised the manuscript. JXL and HZS revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Hui-Zeng Sun.

Ethics declarations

Ethics approval and consent to participate

All samples were collected with the approval of the Animal Care Committee at Zhejiang University (Hangzhou, China). Consent to participate: not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

The changes in the expression (CPM) of fibroblasts marker genes (COL3A1, PAM, and MFAP5) between the newborn and adult rumen tissues from the bulk-seq experiment. Fig. S2. Comparison of epithelial cell landscapes of human stomach and cattle rumen. (A) The UMAP maps of the human stomach single-cell data, cells are colored by cell types. (B) The UMAP maps representing the expression of representative marker genes among epithelial cell types of human stomach. (C) Similarity of epithelial cell types between human stomach and cattle rumen. AUROC scores were used to measure the similarity of cell types: red, high correlation; blue and yellow, low correlation. The AUROC scores in the diagonal are meaningless according to the scoring system and are shown as blanks. Fig. S3. Comparison of epithelial cell landscapes of human skin and cattle rumen. (A) The UMAP maps of the human skin single-cell data, cells are colored by cell types. (B) The UMAP maps representing the expression of representative marker genes among epithelial cell types of human skin. (C) Similarity of epithelial cell types between human skin and cattle rumen. AUROC scores were used to measure the similarity of cell types: red, high correlation; blue and yellow, low correlation. The AUROC scores in the diagonal are meaningless according to the scoring system and are shown as blanks. Fig. S4. The rarefaction of the rumen bacteria based on the 16S rRNA in newborn and adult dairy cattle. NB: newborn; AD: adult. Fig. S5. The relative abundance of the bacterial taxa at the genus level that were significantly enriched in newborn calves and adult cows.

Additional file 2: Table S1.

The information of the number of cells for each cell type. Table S2. Marker genes for each cell cluster for immune cells. Table S3. Marker genes for each cell cluster for non-immune cells. Table S4. Marker genes for each cell cluster of human stomach. Table S5. Marker genes for each cell cluster of human skin. Table S6. Differentially expressed genes and their average read counts in each cell types (adult vs. newborn groups). Table S7. The differential metabolites between newborn and adult groups (adult vs newborn). Table S8. The microbiota-dependent 38 metabolites with high co-occurrence probabilities (the inferred conditional probabilities >1) with Desulfovibrio genus. Table S9. Genes in vitamin B6 binding gene set.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wu, JJ., Zhu, S., Tang, YF. et al. Microbiota-host crosstalk in the newborn and adult rumen at single-cell resolution. BMC Biol 20, 280 (2022). https://doi.org/10.1186/s12915-022-01490-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12915-022-01490-1

Keywords

  • Rumen
  • Host single-cell transcriptome
  • Epithelial microbiota
  • Epithelial metabolome
  • Microbiota-host crosstalk
  • Desulfovibrio
  • Pyridoxal