Age- and Microbiota-Dependent Cell Stemness Plasticity Revealed by Cattle Cell Landscape

Newborn ruminants are considered functionally monogastric animals. The poor understanding of cellular differences between newborn and mature ruminants prevents the improvement of health and performance of domestic ruminants. Here, we performed the single-cell RNA sequencing on the rumen, reticulum, omasum, abomasum, duodenum, jejunum, ileum, cecum, colon, rectum, liver, salivary gland, and mammary gland from newborn and adult cattle. A comprehensive single-cell transcriptomic atlas covering 235,941 high-quality single cells and 78 cell types was deciphered. A Cattle Cell Landscape database (http://cattlecelllandscape.zju.edu.cn) was established to elaborately display the data and facilitate effective annotation of cattle cell types and subtypes for the broad research community. By measuring stemness states of epithelial cells in each tissue type, we revealed that the epithelial cells from newborn forestomach (rumen, reticulum, and omasum) were more transcriptionally indistinct and stochastic compared with the adult stage, which was in contrast to those of abomasum and intestinal tissues. The rapid forestomach development during the early life of calves was driven by epithelial progenitor-like cells with high DNA repair activities and methylation. Moreover, in the forestomach tissues of newborn calves, the Megasphaera genus was involved in regulating the transcriptional plasticity of the epithelial progenitor-like cells by DNA methylation regulation. A novel cell type, the STOML3+ cell, was found to be newborn-specific. It apparently plays a crucial role in stemness maintenance of its own and cholangiocytes in the hepatic microenvironment. Our results reveal that the age- and microbiota-dependent cell stemness plasticity drives the postnatal functional maturity of ruminants.


Introduction
During the Neolithic, ruminants were domesticated to provide meat and milk for humans [1]. They have highly organized digestive systems including forestomach (rumen, reticulum, and omasum), abomasum, small intestines, large intestines, and accessory glands [2,3]. Nutrients absorbed by the gastro intestinal tract (GIT) are delivered to the liver for catalytic metabolic processes to enable cell growth and production. The newborn (NB) ruminant GIT is functionally equivalent to that of monogastric animals and displays marked changes in morphology, structure, and function before reaching full ru minating capacity [4]. For example, the rumen expands grad ually after birth to more than 70% of the GIT volume [5,6]. The maturation of the ruminating system starts after birth and is orchestrated by diverse cell types across multiple tissues. During this process, the digestion undergoes a shift from pri marily intestinally absorbed colostrum and milkderived nutrients to shortchain fatty acids, amino acids, and other compounds from feed and ruminal microbial sources; the metabolic function of the liver changes from being primarily responsible for ketogenesis and glycolysis to gluconeogenesis using propionate as a substrate [4].
The molecular mechanisms behind the postnatal ruminat ing maturation is considered as ontogenic [7], with the mech anisms yet not being completely understood. However, there are recent efforts to link the host transcriptome with microbi ota using lamb or calf models [7,8]. The bovine GIT tissues undergo rapid cell proliferation and differentiation in the post natal period [9], and organogenesis is dependent on cell stem ness of individual cells [10,11]. There is little information available on the unique cellular differentiation markers of bovine cell lineages [12,13], which hinders the determination of stemness states of individual cells and prevents a full under standing of the molecular mechanisms relating to postnatal maturation of the ruminating system. The advances on single cell RNA sequencing (scRNAseq) allow the disentangling of heterogeneity in cell composition and celltypespecific dynamic changes in time of organ development [14,15]. Novel algorithms for quantitatively measuring cellular stemness states at the singlecell level from scRNAseq data have emerged [11] and were successfully applied in a recent study on human organ development by assessing singlecell stemness [16]. In our previous study, we provided novel insights into cell com position and potential cell functions at singlecell resolution in adult (AD) dairy cattle [17]. However, alterations in cell type composition and stemness at multiple temporal and spa tial scales have not been investigated in sufficient detail yet.
The postnatal functional maturity of ruminants is mainly due to the changes in the composition and function of various tissue cell types and microbiota in the digestive system. To investigate the biological mechanisms of postnatal functional maturity in ruminants, we firstly used scRNAseq to build a comprehensive singlecell compendium across 13 bovine tis sue types from the NB and the AD stages and delineated the multitissue patterns of cell stemness plasticity. Next, we inte grated the epithelial microbiome, metabolome, and singlecell transcriptome to investigate the role of microbiota in the cell stemness plasticity of the forestomach. Finally, the new cross species comparison based on singlecell transcriptomic data will extend our knowledge on cattlespecific novel cell types contributing to the maintenance of cell stemness in postnatal maturation of the ruminating system.

Construction of a cattle cell atlas
To delineate the agedependent cell types and transcriptional profiles of various tissues of highly organized digestive systems in cattle, we used 10x Genomics based scRNAseq to establish a profile of forestomach (rumen, reticulum, and omasum), abomasum, small intestines (duodenum, jejunum, and ileum), large intestines (cecum, colon, and rectum), liver, salivary gland, and mammary gland tissues of both NB and AD cattle ( Fig. 1A and Tables S1 and S2). After removing lowquality single cells and doublets (see Materials and Methods), 235,941 highquality cells were retained: 49,689 cells from rumen, 12,844 from reticulum, 20,195 from omasum, 12,637 from abomasum, 16,222 from duodenum, 22,027 from jejunum, 16,001 from ileum, 13,987 from cecum, 12,383 from colon, 11,848 from rectum, 20,066 from liver, 4,894 from salivary gland, and 23,148 from mammary gland (Fig. 1B).
We combined all the 235,941 cells into the cell clustering analysis across tissues with batch effect correction (see Mate rials and Methods) and identified 78 cell clusters (Fig. 1C). The annotation of cell clusters was performed by both data driven analysis and manual interpretation based on the highly expressed marker genes. Detailed information on the descrip tions of cell types and their highly expressed marker genes are available in Table S3 and Fig. S1. A uniform manifold approx imation and projection (UMAP) of the singlecell dataset ( Fig.  1D) revealed that cells were organized primarily by tissue type rather than age. Therefore, to further systematically investigate the dynamics of cell type composition between NB and AD, we performed cell clustering analysis of each tissue type separately (see Materials and Methods). We next explored the cell types that were agedependentpresent in each tissue type. In a single tissue type, any cell type with more than 20 cells in AD cattle and less than 5 cells in NB cattle was considered as ADdependentpresent cell type and otherwise was consid ered as NBdependentpresent cell type. By comparing the numbers of each epithelial or immune cell types between NB and AD stages in each tissue, we discovered the agedependent present cell types (Fig. 1E). For example, MKI67 + Th17 cells were enriched in the AD reticulum and omasum tissues; the chief cell_3 was enriched in the AD abomasum, whereas the ALOX5 + epithelial cell and the chief cell_4 were enriched in the NB abomasum; the BPIFA2C + enterocyte was enriched in the AD rectum, whereas the mast cell was enriched in the NB rectum.
Furthermore, we constructed the Cattle Cell Landscape (CCL) database ( Fig. 2; http://cattlecelllandscape.zju.edu.cn), a comprehensive web resource consisting of 30 singlecell data sets across 13 tissue types. The results of cell clustering, celltype annotation, and marker gene identification in single or multiple singlecell datasets can be available for public access and inter active exploration by using the userfriendly interface in the CCL. The CCL offers users 2 modules, the "Landscape Datasets" module and the "Marker Genes Search" module, to explore the singlecell datasets. In addition to the visualization of the single tissuetype dataset, a comparative analysis of cell types between NB and AD groups is also offered in the CCL. Our scRNAseq analysis revealed special marker genes for distinguishment of different cell types and thus presents an accurate and compre hensive resource of marker genes for cell types across various cattle tissues.

Multitissue representation of cell stemness plasticity paradigm
Previous studies reported that the physiological changes be tween NB and AD cattle were associated with the GIT epi thelium [18,19], and the organogenesis is dependent on cell stemness of individual cells [10,11]. However, little is known about how epithelial cell populations and stemness status of GIT and other key metabolic tissues are changed between the NB and AD stages. Therefore, we focused on the epithelial cellular composition of each organ separately ( Fig. 3A to C and Figs. S2 to S5). Firstly, we compared the epithelial cell differentiation states between NB and AD cattle by calculating the singlecell entropy using the SLICE algorithm [11]. The results showed that there was no significant difference in the entropy of epithelial cells of the jejunum, ileum, liver, salivary gland, and mammary gland tissues between the 2 groups ( Fig.  3D). Interestingly, when compared with the AD group, NB epithelial cells from the abomasum and some gut tissues (duo denum, cecum, colon, and rectum) showed significantly lower entropy (P < 0.0001), while the forestomach (rumen, reticu lum, and omasum) epithelial cells had significantly higher (P < 1.0 × 10 −200 ) entropy in the NB groups (Fig. 3D). Higher entropy is associated with larger cell stemness and functional uncer tainty [11]. These results suggest that the high differentiation potential of forestomach epithelial cells at NB stage is likely the origin of the drastic changes in organ development and functional transition from preruminating to ruminating.
The forestomach tissues undergo rapid growth and devel opment after birth. A previous study found that the most significant transcriptional differences between NB and AD ruminants were observed in the forestomach tissues by explor ing transcriptomes of the entire GIT tissues during the transi tion from preruminant to ruminant [19]. However, which cell types drive the differences in transcriptional plasticity between NB and AD forestomach tissues remains unknown. Therefore, we next examined forestomach tissues to further explore cell stemness at the epithelial cellsubtype level. By comparing the singlecell entropy of the epithelial cell subtypes between   NB and AD groups in the rumen (Fig. 4A), reticulum (Fig. 4B), and omasum (Fig. 4C), we found that almost all the epithelial cell subtypes in the NB group had significantly higher entropy than in the AD group in these 3 tissues. Remarkably, the mitotic cells (MCs; highly expressing MKI67 and KRT14) had the high est entropy among all the epithelial cell subtypes in these 3 tissues. The gene set enrichment analysis (GSEA) also showed that MC had a high potential for cell division ( Fig. 4D and E). The random forest classification revealed that the MC in the reticulum and omasum had similar identities (98.5% and 99.5%, respectively) to those in the rumen ( Fig. 4F and G). By assessing the similarities of epithelial cell subtypes in the fore stomach tissues using the MetaNeighbour analysis [20], we also observed that the gene expression patterns of MC clusters were conserved among the rumen, reticulum, and omasum at both NB and AD stages (area under the receiver operating char acteristics [AUROC] score ≥ 0.95; Fig. S6). In addition to MC clusters, only 2 and 3 cell type clusters were conserved across the rumen, reticulum, and omasum at the NB and AD stages, respectively (AUROC score ≥ 0.95; Fig. S6). In other words, the similarities of most epithelial cell subtypes among the rumen, reticulum, and omasum were not conserved at both NB and AD stages (AUROC score < 0.95; Fig. S6) although the forestomach tissues are composed of the similar stratified squa mous epithelium. Herein, to further unveil the molecular events associated with differentiation states within the MC, the differentially expressed genes (DEGs) were identified in NB and AD cattle. The upregulated DEGs of the MCs in the reticulum and oma sum were enriched in the epithelial cell differentiation and cell proliferation, which is consistent with our previous findings in the rumen tissue [21]. Thus, the MC in forestomach tissue may possess different differentiation patterns in the early post natal period compared with adulthood. To confirm this, we then      Finally, we explored modules of coregulated genes that were differentially expressed in the different paths through the tra jectories. We identified 51 (NB_M1M51) and 56 gene mod ules (AD_M1M56) in the differentiation trajectories of NB and AD groups, respectively ( Fig. 5C and D). The NB_M12 and NB_M23 were enriched in the NB MC (Fig. 5C), and the AD_M11 was enriched in the AD MC (Fig. 5D). Interestingly, compared with each other, AD_M11, NB_M12, and NB_M23 had many unique genes (Fig. 5E). The unique genes of NB_ M12 and NB_M23 were enriched in DNA repair and methyl transferase complex, respectively (Fig. 5F). Gene set score analysis also showed that the MC in the NB group had high DNA repair (Fig. 5G). The higher activities of DNA repair and epigenetic modification may help to explain the high tran scriptional plasticity observed in NB MC.

Microbiota-dependent transcriptional plasticity in the forestomach tissues of NB cattle
In addition to the regulation of the host itself, early develop ment of the rumen was also triggered by commensal microbes [7]. A recent study has suggested that during early life devel opment, the host epigenetic modification (DNA methylation) by symbiotic microbiota represents a fundamental and broad level of regulation [22]. As outlined above, we found that the MC transcriptional plasticity was associated with the activity of methyltransferase complex in the forestomach tissues of NB cattle. Therefore, we next aimed at exploring mucosal microbiota involved in DNA methylation of host cells. Firstly, we accomplished 16S ribosomal RNA (rRNA) gene sequenc ing to analyze the mucosal bacteria in the forestomach tissues at NB and AD stages. We found that the mucosal bacterial profiles of NB and AD were distantly clustered by a principal co ordinates analysis plot from the rumen, reticulum, and omasum tissues (Fig. 6A to C). At the genus level, we found that 5, 5, and 8 bacteria genera showed higher relative abun dances in the rumen, reticulum, and omasum of NB calves, respectively (Fig. S7), of which 3 (Megasphaera, Streptococcus, and Enterococcus) were shared by the rumen, reticulum, and omasum of NB calves.
Streptococcus includes facultatively anaerobic bacteria that could participate in the creation of a reduced environment for anaerobic microbiota [23,24], which might facilitate early col onization of anaerobic microbiota in the forestomach tissues of calves after birth. Emerging evidence has suggested that Megasphaera was associated with the development of calf rumen epithelium [7,8]. Thus, we next performed a metabo lomics analysis in NB and AD forestomach mucosal tissues (see Materials and Methods), and identified 718, 722, and 722 metabolites in the rumen, reticulum, and omasum, respec tively (Tables S4 and S5). Then, we used a microbemetabolite vectors (mmvec) neural network analysis [25] to identify the Megasphaera-host interaction by calculating the cooccurrence probabilities between the Megasphaera genus and metab olites in the forestomach mucosal tissues (see Materials and Methods). We found that 136, 129, and 138 metabolites had high cooccurrence probabilities (the inferred conditional probabilities > 1) with the Megasphaera genus in the rumen, reticulum, and omasum, respectively (Table S6). Only 5 me tabolites, including the dimethyl fumarate, carnitine C8:0, pro pionylcarnitine, NAlphaacetyllysine, and 5methylcytosine, showed high cooccurrence probabilities with the Megasphaera genus and were shared in the rumen, reticulum, and omasum ( Fig. 6D and E). Three of them (propionylcarnitine, NAlpha acetyllysine, and 5methylcytosine) are microbiotadependent metabolites according to a reference library of microbiota dependent metabolites [26]. Finally, we performed a gene set score analysis to calculate the molecular signature scores of methyltransferase complex gene sets of MCs from NB and AD in the forestomach tissues. The results showed that MCs in the NB group had high methyltransferase complex scores com pared with the AD group in the forestomach tissues (Fig. 6F). Thus, the Megasphaera genus may regulate the transcriptional plasticity of forestomach epithelial cells by epigenetic regula tion in NB calves. Further studies to reveal the underlying molecular mechanisms are required.

Age-and species-specific novel cell type enables cell stemness maintenance in NB cattle
In addition to the forestomach tissues, the function of the liver also undergoes dramatic changes during postnatal functional maturation in ruminants. The physiological changes in the liver of NB and AD cattle could be observed and appreciated obviously [4]; however, the driving factors remain poorly understood. After analyzing liver tissues from NB and AD, we found a previously undescribed cell type, the STOML3 + cell, in the liver of NB cattle ( Fig. 7A and B). When the similar ities between liver cell types were compared using the hier archical cluster analysis, the STOML3 + cell clustered much closer to cholangiocytes, hepatocytes, and lymphatic endothe lial cells than to other cell types (Fig. 7A). This cell type has many highly expressed marker genes, such as STOML3, HOXA9, MFAP2, INKA1, and TGFB2 (Fig. 7C), suggesting that it may have a cell stemness maintenance function. For example, HOXA9 is involved in regulating cell stemness [27]. Similarly, INKA1 is a promoter of stem cell latency [28]. By constructing ligand-receptor interaction maps using the CellChat [29], we found that the STOML3 + cell predominantly interacts with itself and a cholangiocyte through KIT signaling (Fig. 7D). KIT signaling has been uncovered to be involved in the differ entiation, proliferation, and migration of many cell types and thereby plays an important role in stem cell maintenance [30]. Collectively, this novel cell type found in NB liver tissue may contribute to the maintenance of cell stemness in the hepatic microenvironment.
Since NB calves are considered functionally monogastric animals and their liver tissues have many functions that were similar to that in monogastric species [4,31], thus, we next explored by crossspecies comparison based on singlecell transcriptomics whether the STOML3 + cell also exists in hu man and mouse liver tissues (see Materials and Methods). We integrated 5 human [32] and 2 mouse liver singlecell datasets [33] with our NB cattle liver dataset and then performed Seurat clustering analysis. The UMAP projection of cells showed a certain extent of cell population variance across cattle, human, and mouse liver tissues (Fig. 7E). In total, we identified 39 distinct clusters (Fig. 7F), and cluster 34 was identified as the STOML3 + cell on the basis of highly expressed marker genes (Fig. 7G). Interestingly, the STOML3 + cell was mainly observed in NB cattle liver tissue and was almost absent in human and mouse liver tissues (Fig. 7H). Although 4 and 5 cells belonging to the STOML3 + cell cluster were observed in the mouse and human liver tissues, respectively, however, they were regarded as mislabeled cells because of the unexpressing marker genes STOML3, MFAP2, INKA1, and, TGFB2. This is a predictable artifact of the annotation scheme because entire clusters in stead of individual cells were annotated in each tissue. Suffi ciently small populations of cells that were algorithmically grouped with a more populous cell type will be misannotated [34]. These results indicate that the STOML3 + cell did not exist in human and mouse liver tissues and could be NB calves specific. Further deeper understanding of the STOML3 + cells is required in the future.   0  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19   20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36 CDCA7L  GMPR  CCR9  TCF4  CIRBP  TGFB2  INKA1  PAN3  GUCY1A1  EML4  ABCG2  GNG11  MFAP2  HOXA9  ID1  FABP5  STOML3   0

Discussion
Constructing agedependent cell atlases for 13 different cattle tissues provided a comprehensive understanding of cell com positions and their dynamic changes. We observed the land marks of singlecell gene expression in different bovine tissues to systematically annotate the cattle cell types and established the CCL website (http://cattlecelllandscape.zju.edu.cn) for free access to substantial celltypespecific gene markers, which may become a valuable reference source in the effective annotation of cell types for the broad research community.
In the life of ruminants, one of the most important physio logical transformations is from the preruminating to the rumi nating stage. Several studies have attempted to identify the underlying molecular mechanisms mediating postnatal matu ration of ruminating calves using bulk RNA sequencing [7,35]; however, these studies only focused on the rumen and offered some opposite findings. As an example, Connor et al. [35] found that the PPARA was a vital gene for the rumen epithelial devel opment during the transition from preruminating to ruminat ing, but in the study of Malmuthuge et al. [7], the PPARA did not exhibit a pattern of temporal expression with calf age. Those differences may be attributed to cell heterogeneity at the tran scriptomic level. Moreover, the complicated biological process from preruminating to ruminating states is likely orchestrated by several tissues. In addition to the rumen, the reticulum, oma sum, abomasum, small and large intestine tissues, and other metabolic tissues undergo wide transformation on the gene ex pression and physiological functions during this process [4,5,18]. In the current study involving the cell stemness of epi thelial cells across 13 cattle tissue types, we observed that the abomasum, duodenum, cecum, colon, and rectum had lower entropy (low cell stemness) in the NB stage compared with the AD stage. The lower entropy represents a skewed distribution of expressed genes in functional gene clusters of individual cells that are associated with welldefined cell fates and functionalities [11]. As such, the abomasum and intestinal tissues (duodenum, cecum, colon, and rectum) may undergo rapid functional spec ificity after birth to be physiologically adapted to colostrum and milk.
Contrary to the results in the abomasum and intestinal tis sues, our results showed that in the early postnatal period, higher levels of cell stemness were measured in the forestom ach (rumen, reticulum, and omasum) epithelial cells. Orga nogenesis of complex organisms is dependent on stemness states of individual cells [10,11]. The high level of cell stemness of the forestomach epithelial cells at the NB stage may be neces sary for the marked morphological and functional changes of these organs after birth. As far as we know, no studies have reported the difference in stemness states between epithelial cells of NB and AD cattle because traditionally the cellular stemness status is measured on the basis of the expression level of known differentiation markers [10,11] or by the cellular morphological features [36], in which the unique cellular differ entiation markers of bovine cell lineages remain under inves tigated [12,13]. On the basis of the newly acquired single cell data, we uncovered that MCs (epithelial progenitorlike cells) in the forestomach tissues (rumen, reticulum, and omasum) had the highest stemness and showed more diverse differen tiation patterns in the NB stage, indicating that the rapid tissue development in the early postnatal life of calves was celltype driven. The higher activities of the epigenetic modification of these cells suggest that the epigenetic status is involved in the process of forestomach postnatal development because each cell has its differentiation patterns that are influenced by its epigenetic status [37]. Taken together, our results offer an agedependent multitissue representation paradigm of cell stemness plasticity, which is associated with the major func tional conversion from the abomasum and gut to the fore stomach that is considered to be the basis of maturation of the ruminating system. Further analyses are required to verify the proposed impacts of cell stemness plasticity observed in our current study and the exact molecular mechanisms.
By integrating microbial metagenomics and host transcrip tomics, a previous study revealed that parts of the host mRNA transcripts were responsive to commensal microbes during the early development of the rumen [7]. In the current study, by integrating epithelial microbiome, metabolome, and singlecell transcriptome, we found that the involvement of the mucosal Megasphaera genus in the proliferation and development of the epithelium may be related to propionylcarnitine, NAlpha acetyllysine, and 5methylcytosine in the rumen, reticulum, and omasum tissues of NB calves. The synthesis of propionyl carnitine requires the participation of lysine and its derivatives, and the propionylcarnitine improves endothelial cell function, in enhancing the blood flow of the blood vessels [38], thereby supporting nutrient supply and utilization. This may closely relate to calf forestomach epithelial tissue growth because the nutrient supply of cells in the forestomach tissues that have not yet matured up to rumination function depends on nutrient transport in blood vessels. DNA methylation plays a vital phys iological role in the regulation of gene expression of cells [39]. The Megasphaera genus had a high cooccurrence probability with 5methylcytosine, and epithelial progenitorlike cells had higher activity of the methyltransferase complex, implying that the Megasphaera genus may participate in regulation of differ entiation and proliferation of epithelial progenitorlike cells in the NB rumen, reticulum, and omasum tissues by epigenetic regulation. In addition, the Megasphaera can convert lactate into butyrate (1 kind of shortchain fatty acid) that is respon sible for modulating the teneleven translocation methylcyto sine dioxygenase enzymatic activity that are involved in DNA methylation by affecting TCA intermediates [22,40]. However, more studies are required to unravel the underlying molecular mechanism.
In the time of the early life stage of calves, the GIT is con tinually communicating with other organs, such as liver [18]. Although the liver of NB ruminants has many functions that were similar to that of monogastric species [31], we found an NB cattlespecific cell type, the STOML3 + cell, which plays a crucial role in stemness maintenance of itself and cholangio cytes. Cholangiocytes are highly dynamic epithelial cells, and their functional physiological modification of hepatocyte derived bile is a functional response of the liver to diets [41]. As such, crosstalks between the STOML3 + cell and the chol angiocyte may contribute to accommodating digestive adapta tions during the transition period. Why STOML3 + cell only exists in cattle liver in the early postnatal period and whether it affects the liver maturation process in response to GIT devel opment warrant further investigations.
There are certain limitations in our current study, which need additional experimental work and new tools. Although a previous study reported that rumen epithelial tissues of calves before and after weaning and corresponding cell clusters were generally similar [42], more tissuetype samples from different time points should be collected and analyzed for an even better understanding of the process. The comprehensive marker genes of the STOML3 + cell were developed in our study. However, we failed to verify the STOML3 + cell by performing the immu nohistochemistry experiments using polyclonal antibodies against humans since useful immunohistochemistry reagents (polyclonal or monoclonal antibodies) with specificity to bovine cells are absent, but we believe that this issue will be solved in near future.
Taken together, we delineated a comprehensive singlecell compendium covering 13 tissue types of NB and AD cattle. We have identified cell heterogeneity for many tissues that have not been well characterized yet and constructed the CCL website to explore in detail the results of singlecell datasets in different bovine tissue types. Moreover, the alterations in cell type com position, cellcell communication, and cell stemness on mul tiple temporal and spatial scales were observed, collectively providing novel and fundamental insights into postnatal rumi nating maturity.

Experimental animals
All the experimental procedures were approved by the Animal Care Committee at Zhejiang University (Hangzhou, China). scRNAseq datasets from AD rumen, reticulum, omasum, abomasum, ileum, rectum, liver, salivary gland, and mammary gland were obtained from our previous study [17]. The AD duodenum, jejunum, cecum, and colon were collected from a single lactating Holstein dairy cow. Three NB calves of similar weight were humanely euthanized in a surgery room to collect the tissue samples. Clinical information on all the animals is available in Table S1.

Single-cell suspension preparation and RNA sequencing
The AD duodenum, jejunum, cecum, and colon samples and the NB reticulum, omasum, abomasum, duodenum, jejunum, ileum, cecum, colon, rectum, liver, and salivary gland tissue samples were freshly collected from the cattle that were humanely euthanized.
After removal of the outer muscle layers, the AD duode num, jejunum, cecum, and colon tissue samples were cut into 1mm pieces and treated with different enzymes for different durations: the duodenum and cecum samples were dissociated using the Lamina Propria Dissociation Kit (Miltenyi Biotec) for 30 min at 37 °C; the jejunum and colon samples were dis sociated using the Multi Tissue Dissociation Kit 2 (Miltenyi Biotec) for 41 min at 37 °C. Then, 10% of fetal bovine serum (FBS) was added to stop the digestion, followed by a step of filtration with the 40μm sieve. Next, after centrifugation at 300 × g for 5 min at 4 °C and after washing with 1× Dulbecco's phosphatebuffered saline (DPBS) containing 2% FBS and 0.2 mM EDTA, the samples were treated with Red Blood Cell Lysis Solution (Miltenyi Biotech) for 2 min. Dissociated cells were washed with 1× DPBS containing 2% FBS and 0.2 mM EDTA, centrifuged at 300 × g for 5 min at 4 °C, and resus pended in 1× DPBS containing 10% FBS and 0.2 mM EDTA.
To prepare the singlecell suspensions of the NB reticulum, omasum, abomasum, small intestines (duodenum, jejunum, and ileum), and large intestines (cecum, colon, and rectum) tissue samples, we referred to the protocols of our previous study [17]. Specifically, we first removed the muscle layers of samples and then chopped the samples into 10 × 0.5 mm 2 pieces. After incubating with 20 mM EDTA for 30 min at 37 °C, samples were rinsed with DPBS and minced into 1mm pieces. Tissue samples were incubated in a 37 °C water bath for 5 min with 0.25% TrypsinEDTA (Gibco) and then placed in ice for 2 min, while the digestion was stopped by adding the prechilled Hank's balanced salt solution (HBSS). Samples were centrifuged at 4 °C for 2 min with 300 × g and removed the supernatant and then washed twice with cold HBSS. Next, for the NB retic ulum, omasum, and abomasum samples, the multiple enzymes (1.5 mg/ml collagenase IV, 1.5 mg/ml collagenase I, 100 U/ml hyaluronidase, 50 U/ml deoxyribonuclease I [DNase I], and 1.5 mg/ml dispase) were added and treated for 30 min at 37 °C. For the NB small intestines (duodenum, jejunum, and ileum) and large intestines (cecum, colon, and rectum) tissue sam ples, the multiple enzymes from the Propria Dissociation Kit (Miltenyi Biotec) were added and treated for 30 min at 37 °C. The digestion was terminated by adding 10% FBS, followed by a step of filtration with the 70 and 30μm SmartStrainer (Miltenyi Biotec). After centrifugation at 300 × g for 5 min at 4 °C, samples were resuspended in 2 ml of HBSS. After centrif ugation at 4 °C for 5 min with 300 × g and washing twice with 1× phosphatebuffered saline (PBS) with 0.04% bovine serum albumin (BSA). Finally, samples were centrifuged at 4 °C for 5 min with 300 × g and resuspended in 1× PBS with 0.04% BSA.
Singlecell suspensions of the liver and salivary gland tis sues of NB calves were prepared according to the protocol of our previous study [17]. For the NB liver tissue, the sample was transferred into cold DPBS, chopped into 1mm pieces, and treated with 1.5 mg/ml collagenase I, 20 U/ml papain, and 50 U/ml DNase I for 20 min at 37 °C. The dissociated cells were then passed through 70 and 30μm SmartStrainer, cen trifuged at 4 °C for 5 min with 300 × g, and resuspended with 2 ml of cold HBSS. Dissociated cells were centrifuged for 5 min at 300 × g, washed twice with 1× PBS with 0.04% BSA at 4 °C, and then resuspended in 1× PBS with 0.04% BSA. For the NB salivary gland tissue, the sample was isolated from the parotid gland was placed into cold 1× PBS, and then the blood and fat was removed. The sample was cut into 0.5mm pieces and incubated with 1.5 mg/ml collagenase II, 1.5 mg/ml dis pase, and 50 U/ml DNase I at 37 °C for 20 min. Dissociated cells were passed through 70 and 30μm SmartStrainers. After centrifugation at 4 °C for 5 min with 300 × g, dissociated cells were resuspended with 300 μl of RPMI 1640 medium (Gibco) and treated with 3 ml of Red Blood Cell Lysis Solution for 3 min. Then, dissociated cells were centrifuged at 4 °C for 5 min with 300 × g and resuspended with 2 ml of RPMI 1640 medium. After centrifugation for 5 min with 300 × g, disso ciated cells were washed twice with 1× PBS with 0.04% BSA at 4 °C and then resuspended with 1× PBS with 0.04% BSA.
The Countess II Automated Cell Counter was used for the assessment of the viability of singlecell suspension via trypan blue. If the cell viability was low, then the MACS Dead Cell Remove Kit (Miltenyi Biotec) was used to remove dead cells according to the manufacturer's recommendations. Dissociated live cells were captured to construct the libraries using the Chromium Single Cell 3' Reagent Kits v3 (10x Genomics) following the manufacturer's recommendations. The Agilent Bioanalyzer High Sensitivity chip was used to check the quality of the libraries, and then libraries were sequenced on the NovaSeq 6000 platform (150bp pairended manner).

scRNA-seq data analysis
We generated 15 new scRNAseq datasets in the present study. The CellRanger (version 3.1.0) was used for sample demulti plexing, barcode processing, and singlecell 3' gene counting. The reads of scRNAseq data were aligned to the ARSUCD1.2 cattle reference genome. Singlecell transcriptomics datasets, generated using 10x Genomics, of AD rumen, reticulum, omasum, abomasum, ileum, rectum, liver, salivary gland, and mammary gland as well as the NB rumen and mammary gland were collected from our a previous study [17] and inhouse lab datasets. A total of 30 scRNAseq datasets representing multiple tissues were analyzed individually. For each dataset, highquality cells with 500 to 4,000 genes detected, unique molecular identifier counting < 50,000, and the ratio of mito chondrial gene < 40% were retained using the Seurat (version 4.0.3) R package [43] for further analysis. The doublets in each dataset were also removed using the DoubletFinder package (version 2.0.3) [44].
To create the crosstissue temporal cell atlas, the aforemen tioned filtered, nonscaled Seurat objects of NB and AD cattle datasets were imported and merged into a Seurat object. The merged object was normalized (NormalizeData function of Seurat) and scaled to regress out the cell cycle signature (com puted using Seurat's CellCycleScoring function and scaled using Seurat's ScaleData, vars.to.regress = c("S.Score", "G2M. Score")). The principal component analysis was performed (RunPCA function of Seurat), and next, the batch effect was corrected using the Harmony R package (version 0.1.0) [45]. The FindClusters function of Seurat with an appropriate res olution was used to accomplish cell clustering. The UMAP embedding algorithm (RunUMAP function) or tdistributed stochastic neighbor embedding algorithm was used to vis ualize cells. Finally, the FindAllMarkers function of Seurat was performed to identify marker genes (|"avg_logFC"| > 0.25 and "p_val_adj" < 0.05).
To further investigate the dynamics between the NB and AD groups, we merged the nonscaled Seurat objects of the NB and AD datasets for each tissue type and performed cluster analysis separately following the same analysis pipeline (but not regressing out the cell cycle signature) to generate age dependent cell atlases for each tissue type. We further manually removed cell clusters (potential doublets) that had overlapping gene profiles of multiple different cell types and excluded the hemocytes. The rumen agedependent cell atlas was collected from our previous study [21].
The singlecell entropy (scEntropy) of individual epithelial cells in NB and AD cattle tissue was quantitatively measured on the basis of scRNAseq datasets using the SLICE pipeline [11] with some modifications. In brief, using Kappa statistics [46], we calculated the functional similarity of genes in the Bos taurus GOTERM_BP_FAT subset that was collected from the DAVID Functional Annotation dataset [47] to generate the genomewide genetogene functional similarity matrix of Bos taurus. Next, we replaced the functional similarity matrix of human with that of Bos taurus, and other analysis steps were performed with default parameters following the SLICE pipeline. The differences between NB and AD cattle were assessed using the analysis of variance, with a P value < 0.0001 as the threshold for significance. *, **, and *** represent P < 0.0001, P < 1.0 × 10 −100 , and P < 1.0 × 10 −200 , respectively.
To assess the similarities of epithelial cell subtypes in the forestomach tissues, the expression matrices (raw counts) of epithelial cell subtypes of the rumen, reticulum, and omasum at NB and AD cattle were merged, respectively. The 2 merged expression matrices were separately used to perform the MetaNeighbour analysis [20] with default parameters. The epi thelial cell subtypes were considered to be conserved if AUROC scores between them were higher than 0.95.
We used the Monocle3 [48] to model differentiation tra jectories of reticulum epithelial cell types in NB and AD cattle with default parameters following the general pipeline (https:// coletrapnelllab.github.io/monocle3/). The CellChat [29] with default parameters was used to perform ligand-receptor inter action analysis to recognize the inferred cellular communication patterns in the NB liver scRNAseq data. The "AddModuleScore" function of Seurat (version 4.0.3) was used to perform gene set scoring analysis to calculate the signature scores on the gene sets in the MCs between NB and AD cattle in the fore stomach tissues. The 2sided Wilcoxon rank sum test was used to evaluate the differences in the signature scores between NB and AD cattle. Detailed information of the gene sets "DNA repair" and "methyltransferase complex" are listed in Table S7.
To explore the STOML3 + cell across different species, we downloaded the dataset from 5 human liver singlecell datasets [32] (https://github.com/BaderLab/HumanLiver) and 2 mouse liver singlecell datasets [33] (GSM3714747 and GSM3714748 in the National Center for Biotechnology Information Gene Expression Omnibus accession GSE129516). We firstly per formed expression matrix preprocessing separately for the 3 species using the Seurat (version 4.0.3) and retained the detected homologous genes of these 3 species in each scRNA seq dataset. Next, we merged the datasets into a Seurat object and performed cell clustering analysis to come up with 39 cell clusters.

CCL website construction
The main CCL website was completed by PHP language and MySQL. The CCL website offers 2 function modules: the "Landscape Datasets" module and the "Marker Genes Search" module. The "Landscape Datasets" module enabled the visu alizations for cell types of 13 different tissue types and allowed investigation the celltype populations' heterogeneity from NB and AD cattle. The "Marker Genes Search" module provided the top 100 highly expressed genes of each cell type across various tissues and also displayed the expression of marker genes presented as a heatmap at the celltype averaged level. More details about the CCL website are available at http:// cattlecelllandscape.zju.edu.cn.

Differential expression analysis
The Wilcoxon ranksum test as implemented in the Seurat's "FindMarkers" function [43] was performed to identify DEGs of MC between NB and AD cattle. Only genes with |LogFC| > 0.5 and adjusted P value < 0.05 were considered as DEGs.

Enrichment analysis
The gene sets that are enriched in MCs of the reticulum and omasum were identified by the GSEA using the "GSEA" func tion of clusterProfiler R package [49]. The "enrichGO" function of clusterProfiler R package was applied to the performance of Gene Ontology (GO) term enrichment analysis based on the dataset "org.Bt.eg.db".

Random forest classifier
The ClassifyCells function that is implemented in the Seurat R Package (version 2.2.0) was used to perform the random forest classification with default parameters. We grouped the subtypes of MC identified from the rumen samples as the training class, which was subsequently applied to the MC from the reticulum and omasum samples.

16S rRNA gene sequencing and analysis
Six animals (3 per group) were selected to collect the mucosa. The mucosal bacterial profiles of the rumen in the NB and AD cattle were collected from our previous study [21]. Total DNA of the mucosal microbiota was extracted from each reticulum and omasum epithelial tissue using the E.Z.N.A.® DNA Kit (Omega BioTek) according to the manufacturer's instructions. The primer pairs 338F (5′ACTCCTACGGGAGGCAGCAG3′) and 806R (5′GGACTACHVGGGTWTCTAAT3′) were used to amplify the hypervariable region (V3V4) of the bacterial 16S rRNA gene. After polymerase chain reaction amplification, all amplicon libraries were sequenced on an Illumina MiSeq platform using the pairedend 2 × 300bp protocol. The paired end reads were merged using the FLASH (version 1.2.11) [50]. The data were processed using the QIIME2 software package [51]. The DADA2 plugin [52] was used to perform further quality control and chimera removal and to produce the ampli con sequence variant (ASV) feature tables. The representative sequences of these ASVs were identified and assigned to the bacteria database of SILVA (Release138, http://www.arbsilva. de) and clustered at 99%. We performed linear discriminant analysis effect size based on the ASVs that were merged to the genus level to identify bacterial taxa that were significantly (linear discriminant analysis > 3.5, P < 0.05) enriched in NB or AD cattle of the reticulum and omasum.
To predict the relationship between the genus Megasphaera and metabolites in the rumen, reticulum, and omasum tissues, the microbe-metabolite vectors (mmvec) neural networks anal ysis was performed [25]; the analysis pipeline using default parameters was performed to determine the conditional prob ability that each metabolite was present given the presence of genus Megasphaera and based on the microbial sequence counts and the metabolite relative concentration.

Metabolomics analysis
The same animals were selected to conduct metabolomics analysis. The mucosal metabolite profiles of the rumen in the NB and AD cattle were collected from our previous study [21]. The epithelial tissues of the reticulum and omasum were homogenized with 1,000 μl of 70% (v/v) icecold methanol/ water and cold steel balls at 30 Hz for 3 min and then were whirled for 1 min without steel balls. After resting for 15 min and centrifugation at 4 °C, 12,000 rpm for 10 min, the super natant of the sample was collected for liquid chromatography tandem mass spectrometry analysis. We performed the electrospray ionization quadruple traptandem mass spectrometry analysis for the supernatant of each sample. The operation parameters of the electrospray ionization source were as follows: ion spray voltage, 5,500 V (positive) and −4,500 V (negative); source temperature, 500 °C; ion source, gas I and gas II; curtain gas were set at 50, 50, and 25 psi, respectively; ion spray volt age, 5,500 V (positive) and −4,500 V (negative); the collision gas was set to high. The instrument tuning and mass calibra tion were accomplished with 10 and 100 μmol/l polypropylene glycol solutions in QQQ and LIT modes, re spectively. The Analyst software (version 1.6.3) was used to analyze the mass spectrometric data.