Single-cell transcriptome and metagenome profiling reveals the genetic basis of rumen functions and convergent developmental patterns in ruminants

Abstract

The rumen undergoes developmental changes during maturation. To characterize this understudied dynamic process, we profiled single-cell transcriptomes of about 308,000 cells from the rumen tissues of sheep and goats at 17 time points. We built comprehensive transcriptome and metagenome atlases from early embryonic to rumination stages, and recapitulated histomorphometric and transcriptional features of the rumen, revealing key transitional signatures associated with the development of ruminal cells, microbiota, and core transcriptional regulatory networks. In addition, we identified and validated potential cross-talk between host cells and microbiomes and revealed their roles in modulating the spatiotemporal expression of key genes in ruminal cells. Cross-species analyses revealed convergent developmental patterns of cellular heterogeneity, gene expression, and cell-cell and microbiome-cell interactions. Finally, we uncovered how the interactions can act upon the symbiotic rumen system to modify the processes of fermentation, fiber digestion, and immune defense. These results significantly enhance understanding of the genetic basis of the unique roles of rumen.

Figure 1.

UMAP plots and the number of cell types. Cell number (n > 3) at the 17 time points (middle) and expression levels of representative marker genes shown as violin plots (right) for each cell type in sheep (A) and goat (B) rumen.

Figure 2.

Immunofluorescence staining and pseudotime trajectory analysis. (A) Immunofluorescence of KRT15 in rumen tissues of sheep at the time points e135 and d21 and of goats at the time points e120 and d7. Scale bars, 60 μm (zoomed-in image). (B) Immunofluorescence of LYVE1 in rumen tissues of sheep and goats at the time points e90 and d90. Scale bars, 60 μm (zoomed-in image). (C) Pseudotime trajectory analysis corresponding to the differentiation of basal cells in sheep and goats. Cells are colored by pseudotime, cell types, and sample stages. (D) Pseudotime trajectory analysis corresponding to the marker gene CA1 in sheep and CA3 in goats. Cells are colored by sample stage.

Figure 3.

scRNA-seq-based gene expression in rumen cell atlases of sheep and goats. (A) Heatmap showing the top 10 expressed genes (row) in each cell (column) in sheep and goat rumen. The right panel shows the significant (P_adj < 0.05) Gene Ontology (GO) terms for the marker genes in each cell type. (B) Heatmaps showing the up-regulated and down-regulated DEGs for each pairwise comparison between different stages (stages E2 vs. E1, E3 vs. E2, P1 vs. E3, P2 vs. P1, and P3 vs. P2) in sheep and goat rumen. The abbreviations of the cell names are shown in Figure 1.

Figure 4.

Changes in core regulatory transcription factors (TFs) during rumen development. (A) Common up-regulated and down-regulated TFs across the six major developmental stages. (B) Heatmap of the area under the curve (AUC) scores of TF motifs estimated per cell of the cell types across the six major developmental stages. The abbreviations of the cell names are shown in Figure 1.

Figure 5.

Changes in ligand–receptor (LR) interactions between different cell types in the comparisons of adjacent stages (stages E2 vs. E1, E3 vs. E2, P1 vs. E3, P2 vs. P1, and P3 vs. P2) in sheep and goats. (A) Functional enrichment analysis showing the GO terms for the common up- or down-regulated LR interactions (P_adj < 0.05). (B) Common interaction pairs and their frequencies (Freq).

Figure 6.

Conserved TFs and cell communications in rumen or stomach tissues across species of sheep (including download data), goats, cattle, monkeys, mice, and humans. (A) Dot plots showing the cell–cell interaction pairs in the indicated cell types. The size of the dots indicates the −log₁₀ P-value, and the colored dots show the mean expression levels. (B) Heatmap of the area under the curve (AUC) scaled scores of TF motifs estimated per cell of the cell types.

Figure 7.

Ruminal microbiota correlations with single-cell transcriptome sequencing in sheep and goats, and transcriptomic analysis of co-cultured rumen cells with Prevotella copri. (A) Correlations of host rumen cell types with the dominant ruminal microbiota. Heatmap shows the Spearman's correlation coefficient between the Eigengene/PC1 value of cell types and the abundances of dominant ruminal microbiota (P < 0.05). (B) Correlations of host rumen cell types with microbial KEGG pathways. Heatmap shows the Spearman's correlation coefficient between the Eigengene/PC1 value of cell types and the abundances of selected microbial KEGG pathways (P < 0.05). The abbreviations of the cell names are shown in Figure 1. (C,D) GO terms (P < 0.05) for up- and down-regulated DEGs (P < 0.05 and |logFC| > 0.75) between the treatment group (co-cultured with P. copri) and control group in primary rumen epithelial cells (C) and primary rumen Fib (D). (E,F) Boxplots showing expression levels of overlap between up-regulated DEGs (treatment vs. control) in primary rumen epithelial cells (E) and primary rumen Fib (F) with up-regulated DEGs (stages P2 vs. P1; P_adj < 0.05 and |logFC| > 0.25), which were significantly (P < 0.05) correlated with P. copri. (****) P < 0.0001, (**) P < 0.01, (*) P < 0.05.

Figure 8.

A summary model of the genetic basis underlying the rumen functions along the development stages.