Scientific Papers

The early life immune dynamics and cellular drivers at single-cell resolution in lamb forestomachs and abomasum | Journal of Animal Science and Biotechnology


Forestomachs share similar gene expression patterns during early developmental stages

To understand the mucosal gene expression patterns of four stomachs of ruminants at early ages, we collected 48 mucosal tissue samples from the rumen, reticulum, omasum, and abomasum, from each of the triplicate lambs at 5, 10, 15, and 25 days of age, and subsequently performed transcriptomic sequencing. After annotation and quality control, a total of 13,174 ± 448 genes (counts per million reads, CPM > 1) expressed in each sample were obtained. The gene expression pattern in abomasum clustered independently from that in forestomachs along the first principal component (PC1), accounting for 36.65% of the total variation (Fig. 1A). Spearman correlation analysis for all pairs of RNA-seq samples (Fig. 1B) further revealed the close similarities in all forestomach samples and the divergence between the abomasum and forestomachs. D 10 and 15 were the separatrices that divided the gene expression pattern in the same tissue into two closer stage blocks (d 5 to 10 and d 15 to 25).

Fig. 1
figure 1

Gene expression patterns of four stomachs during early development. A Principal component analysis (PCA) of 48 samples across four ages and four stomachs based on the gene expression levels. The samples were clustered by tissue (left) and tissue plus age stage (right), respectively. B Heatmap showing the spearman correlation for 48 samples based on the gene expression level. C Bar plot showing the differentially expressed genes (DEGs) between adjacent age groups. D Heatmap showing the relative expression of all DEGs between adjacent age groups, rows (each representing a DEG) and columns (each representing one sample) are unsupervised clustered

Next, we determined the DEGs within the same tissue between any two adjacent ages (d 10 vs. d 5, d 15 vs. d 10, and d 25 vs. d 15). Stages from d 10 to 15 showed the highest number of DEGs in all tissues, followed by stages from d 5 to 10 and d 15 to 25 (Fig. 1C). Unsupervised hierarchical cluster analysis [33] across the ages of Hu sheep (Fig. 1D) showed that the DEGs could be clustered by both age and tissue type. The forestomachs and abomasum were classified into two clusters, and ages from d 5 to 10 and d 15 to 25 were classified into two clusters in the forestomachs and abomasum, respectively. Together, these selected DEGs provided a molecular signature for the developmental status during early developmental stages.

Immune and defense responses in the mucosa of four stomachs are quickly developing with age during early life

To characterize the dynamic changes in gene expression, we clustered all expression patterns using the WGCNA method. Except for the grey module that contained unclustered genes, we identified 10 main gene transcriptional modules, with the turquoise module (Mturquoise) containing the highest number of genes, followed by Mblue, Mbrown, and Myellow (Fig. 2A). We calculated the eigengene of each module and assessed module correlations with age using the module-eigengene. Five modules (Mred, Mblack, Mpink, Mpurple, and Mmagenta) were positively correlated with age, whereas three modules (Mgreen, Myellow, and Mturquoise) were negatively correlated with age (P < 0.05, Fig. 2B). GO and KEGG pathway enrichment showed that among the eight modules, no pathways were enriched in Mblack, Mpurple, Mgreen, or Myellow. Therefore, we focused on the four remaining age-associated modules, where eigengene expression of Mred and Mpink showed increased dynamics in each stomach (Fig. 2C–D), Mmagenta showed increased dynamics only in the abomasum (Fig. 2E), and Mturquoise showed a dramatic drop synchronously from d 10 to 15 in each stomach (Fig. 2F). Mred was associated with cell adhesion and its regulation, as well as immune function, such as “Adaptive immune response”, “Positive regulation of cell–cell adhesion”, “Positive regulation of T cell activation” (Fig. 2G), and “Th1 and Th2 cell differentiation”, “Th17 cell differentiation”, “Hematopoietic cell lineage”, “Cell adhesion molecules”, and “Cell adhesion molecules” (Additional file 2: Fig. S1A); Mpink was involved in defense response, and inflammatory or stimulus–response related pathways, respectively, such as “Defense response to virus” and “Defense response to symbiont” (Fig. 2H), and “RIG-I-like receptor signaling pathway”, “Toll-like receptor signaling pathway”, and “TNF signaling pathway” (Additional file 2: Fig. S1B); Mmagenta was associated with cell cycle, division, and their regulation, such as “Cell cycle”, “Regulation of mitotic cell cycle”, and “Cell division” (Fig. 2I), and “Homologous recombination” (Additional file 2: Fig. S1C); Mturquoise was enriched for peptide biosynthetic and metabolic process, and oxidative phosphorylation, such as “Peptide biosynthetic process”, “Translation”, and “Amide biosynthetic process” (Fig. 2J) and “Ribosome” and “Oxidative phosphorylation” (Additional file 2: Fig. S1D). The key genes driving the function of Mred included CORO1A, CD3E, CD4, TBC1D10C, LAPTM5, IL2RG, MAP4K1, SEPT1, TNFRSF18, and PTPRC (Fig. 2K); of Mpink, included ZBP1, IFI44L, HERC6, DDX58, ZNFX1, XAF1, IFIT1, IFIT3, LOC101103965, and MX1 (Fig. 2L); of Mmagenta, included KIF11, PLK4, ECT2, TOP2A, TTK, NDC80, DLGAP5, MIS18BP1, HMMR, and ATAD5 (Fig. 2M); and of Mturquoise, included ribosomal protein large subunit family genes (including RPL23, RPL7, RPL27A, RPL26, and RPL34), and GTPBP6, TMEM185A, LOC780463, RASSF1, FBXW5, MRPS2, WDR5, and RHOT2 (Fig. 2N).

Fig. 2
figure 2

Gene co-expression network and their function in four stomachs during early development. A WGCNA cluster dendrogram and module assignment. Modules corresponding to branches are labeled with colors indicated by the color bands underneath the tree. The bar plot under the module assignment showing the gene number of each module. B Pearson correlation between modules and ages calculated based on the eigengene expression. C–F Boxplot showing the average eigengene expression of each tissue in Mred (C), Mpink (D), Mmagenta (E), and Mturquoise (F). G–J Top GO pathway enriched in Mred (G), Mpink (H), Mmagenta (I), and Mturquoise (J). Only GO terms for “Biological Process” were selected. K–N Co-expression network of Mred (K), Mpink (L), Mmagenta (M), and Mturquoise (N) constructed based on the WGCNA connection. Each module co-expression network were color-coded in module color. Each node represents one gene, the hub genes (top 10 in Mred, Mpink, and Mmagenta or top 30 in Mturquoise) were at the outermost of each networks

To validate the functional development of the four stomachs, we compared the functional differences using GSEA based on all expressed genes (CPM > 1). We identified 50 KEGG pathways exhibiting significant differential enrichment between adjacent ages and at least in one age stage in one stomach (Fig. 3A). Consistent with the functional enrichment and key genes shown in turquoise, the ribosome pathway was downregulated from d 10 to 15 in all stomachs (Fig. 3A and D), indicating downregulated protein synthesis [34]. However, most of the upregulated pathways were associated with immune function, such as cytokine-cytokine receptor interactions and Th17 cell differentiation, from d 5 to 10 in the rumen (Fig. 3B–C). To better understand dynamic immune function, we constructed an immune-associated gene set containing 34 pathways based on the KEGG database, and compared the activities of these pathways using GSVA. We observed a general increasing trend in the 34 immune-associated pathways from d 5 to 25, in all four stomachs (Fig. 3E). However, these pathways showed a delayed increase (d 15 or 25) in the reticulum and omasum compared to that in the rumen and abomasum (d 5 or 10). Taken together, the results from both GSEA and GSVA suggested a marked functional focus on the immunity and defense of the four stomachs during early sheep developmental stages.

Fig. 3
figure 3

Function dynamics of four stomachs during early development. A Dot plots showing the activities of KEGG pathways enriched by get set enrichment analysis (GSEA) across the adjacent ages (adjusted P value < 0.05). B–D Examples of specific KEGG pathways of A including “Cytokine-cytokine receptor interaction” in rumen between d 10 vs. d 5 (B), “Th17 cell differentiation” in rumen between d 10 vs. d 5 (C), and “Ribosome” in omasum between d 15 vs. d 10 (D). E Heatmap showing immune associated pathway activities of each sample scored by gene set variation analysis (GSVA) across four ages

Immune cells and endothelial cells play essential roles in immune and defense function during early developmental stages

Assessing the immune functions of various cell types is essential to fully elucidate the mechanisms underlying immune function development. In this study, we performed a single-cell sequencing analysis of cells dissociated from the entire mucosa of the rumen, reticulum, omasum, and abomasum of three Hu lambs at 25 days of age. Among these, 5,820, 6,096, 6,420, and 4,433 cells originated from the rumen, reticulum, omasum, and abomasum, respectively (Fig. 4A). We classified these cells into groups of cell types using graph-based clustering, based on informative principal components. Based on the expression of the marker genes (Fig. 4B, Additional file 1: Table S1), we classified these cells into known cell lineages: basal cells, granule cells, spinous cells, proliferative cells, pit/gland mucus cells, chief cells, parietal cells, fibroblasts, endothelial cells, smooth muscle cells, monocytes, macrophages, and T cells (Fig. 4B). The proportion of each cell lineage varied significantly among the four tissue types (Fig. 4C).

Fig. 4
figure 4

Single cell view and signature of four stomachs. A tSNE analysis of 22,589 cells from rumen, reticulum, omasum, and abomasum at d 25, with four tissues labeled in left side, and 12 major cell types labeled in right side with different colors. B Dot plot showing the expression of marker genes for the cell types. C The proportion of each cell type among four different tissues. D Dot plot showing the expression of hub genes of Mred, Mpink, Mmagenta, and Mturquoise in 13 cell types; the hub genes of each module were identified through WGCNA analysis (see Fig. 2E). E Average expression of marker genes for monocytes/macrophages, T cells, proliferative cells, and endothelial cells of four tissues across four ages in the transcriptomic data, respectively. Data across four ages were analyzed using one-way analysis of variance with Tukey’s multiple-comparisons test (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, ns: no significant)

To further identify the major cell types dominating the maturation of immune and defense functions, we investigated the expression of the hub genes identified in Mred, Mpink, Mmagenta, and Mturquoise for each cell type (Fig. 4D). All of the Mred hub genes were highly expressed in T cells; among them, CORO1A, LAPTM5, IL2RG, MAP4K1, and PTPRC were highly expressed in both T cells and monocytes or macrophages. Eight Mpink hub genes (IFI44L, DDX58, ZNFX1, XAF1, IFIT1, IFIT3, LOC101103965, and MX1) were highly expressed in endothelial cells. Consistent with the module enrichment analysis results, all magenta hub genes were highly expressed in proliferative cells. However, the expression of the turquoise hub genes did not significantly differ among the different cell types. These results suggest that T cells, monocytes, or macrophages dominate Mred; endothelial cells and proliferative cells dominate Mpink and Mmagenta, respectively; and Mturquoise is dominated by multiple cell types that work together.

We assessed the average relative expression of marker genes (the marker genes were consistent with Fig. 4B) for monocytes/macrophages, T cells, and endothelial cells, in the transcriptomic data across the four time points. We observed an age-dependent increase in the expression of monocytes/macrophages and T cells marker genes in the four stomachs (Fig. 4E). There was an increase in only the proliferative cells from d 5 to 10 in forestomachs, and a continuous increase in proliferative cells in the abomasum, indicating the ratio dynamics of the above cell types from d 5 to 25. However, the expression of endothelial cell marker genes did not change over time in the rumen and reticulum, and it decreased from d 10 to 25 over time in omasum and abomasum. This phenomenon combined with the dynamic of Mpink hub gene expression indicated that the expression of Mpink hub gene of endothelial cells rather than the cell ratio increase with age.

The integration of single-cell results and functional enrichment analyses through GSVA and GSEA confirmed that the main events of development were age-dependent increases in immune cells and their immune functions, as well as in endothelial cells and their defense-related functions in the early developmental stages of Hu sheep.

Widespread immunity activation occurs in different stomachs

In addition to classical immune cells, non-immune cells play an important role in the immune function of different tissues [35]. To determine whether non-immune cells in four stomach tissues are involved in immune functions, we systematically investigated the global involvement of all cell types in individual tissue immunity. Based on the known cell markers (Additional file 2: Fig. S2A–D; Additional file 1: Table S2–S5), we defined 19, 19, 16, and 18 clusters of rumen, reticulum, omasum, and abomasum cells, respectively (Fig. 5A). The forestomachs share similar cell types including basal cells (BC), spinous cells (SC), granule cells (GC), proliferative cells (PC), fibroblasts (FC), endothelial cells (EndoC), T cells (TC) and monocytes/macrophages (MC). In addition, a type of smooth muscle cells (SMC) in reticulum and undefined epithelial cells (UEC) in omasum were identified. In addition to PC, FC, EndoC, TC and MC, the diverse cell types including parietal cells (ParC), chief cells (CC), pit mucus cells (PMC), gland mucus cells (GMC) and enteroendocrine cells (EnteC) were identified in abomasum.

Fig. 5
figure 5

Widespread activation of immune function and their associated gene expression in different cell types in four stomachs. A t-SNE maps of the cells in rumen, reticulum, omasum, and abomasum (left to right). B Immune-associated pathway activities of each cell type in rumen, reticulum, omasum, and abomasum (left to right). C Dot plot showing the expression patterns of genes associated with the defense or response to microorganism in different cells in rumen (upper left), reticulum (upper right), omasum (lower left), and abomasum (lower right). The genes and their function were obtained from GO database, only genes that expressed in at least one cell type (Percent expressed > 0) were kept. The names of the GO terms were used as the summary of the function of these genes ()

We performed a cross-tissue comparison of all cell types for immune- and defense-related pathway enrichment using GSVA. In addition to immune cells (TC and MC) in all four tissues, some non-immune cell types had high immune activities, particularly EndoC, in the four stomachs. In addition, GC3 in the rumen, FC2 and SC1 in the reticulum, UEC and GC in the omasum, and PMC2 and FC in the abomasum exhibited relatively high immune activity (Fig. 5B). These results suggest that non-immune cells undergo a layer of cellular regulation for tissue-specific immunity. Among these immune-related pathways, the toll-like receptor and NOD-like signaling receptor signaling pathways exhibit distinct activities among different cell types due to their important roles in both innate and adaptive immune activation [36, 37]. Although immune cells are central to pattern recognition receptor (PRR) signaling responses, many other cell types, including epithelial, stromal, and endothelial cells, express PRRs [36], which recognize a broad range of microbial-associated molecular pattern ligands and trigger immune activation. Therefore, we further explored the expression of PRRs, including that of toll-like receptors (TLRs), NOD-like receptors (NLRs), and RIG-I-like receptors (RLRs), in all cell types across the four stomach tissues (Additional file 2: Fig. S3A–D). Notably, there is a distinct group of non-immune cells, including GC3 in the rumen, EndoC1 in the reticulum, FC in the omasum, and EndoC in the abomasum, that show high expression of PRR genes. We also investigated the expression of G protein-coupled receptors (GPCRs) that mediate most of the physiological responses to various signaling molecules including microbial metabolites, thus involved in regulating gastrointestinal mucosal immunity and maintaining intestinal barrier (Additional file 2: Fig. S4A–D). Here we focused on fatty acid receptors because they are most studied and abundant group of microbial metabolites, we also focused on the members of GPCR families that have been validated to be activated by microbiome cultures in a recent high-throughput screening study [38]. The GPCR families shared similar expression patterns between the mucosa of forestomachs, compares to those in abomasum. For example, the free fatty acid receptor 2 (FFAR2, also known as GPR43) and free fatty acid receptor 3 (FFAR3, also known as GPR41) that sense short-chain fatty acids (SCFAs), were lowly expressed in all cell types, whereas the FFAR4 (also known as GPR120) that senses medium- and long-chain unsaturated fatty acids, was highly expressed in the GC and SC and some subtypes of BC in forestomachs. In addition, fatty acid binding protein 4 (FABP4) and 5 (FABP5) were highly expressed in most of PC, BC, SC and GC subtypes in forestomachs. However, FABP3 was highly expressed in ParC in abomasum. Although succinate receptor (SUCNR1), adrenergic receptors (ADRA family), cholinergic receptors (CHRM family), dopamine receptors (DRD family), histamine (HRH family) and 5-hydroxytryptamine receptors (HTR family) can be activated by multiple microbial strains culture [38], most genes coding for these families were lowly expressed. For example, SUCNR1 that can be activated by various microbial culture including the strains belong to Bacteroides, Prevotalla, Escherichia coli, were not nearly expressed in the all cell subtypes in forestomachs and abomasum. We also noted that F2RL1, HTRA1, HTRA2 were highly expressed in the most of PC, BC, SC, GC, FC and EndoC subtypes in forestomachs.

In addition to the differences in the receptors that recognize microorganisms and their metabolites, the key responses of different cells to various microorganisms remain undefined. Therefore, we analyzed the expression of genes related to host defense and responses to bacteria, fungi, protozoans, and viruses, based on the GO database (Fig. 5C). Genes related to defense and response to bacteria were highly expressed in epithelial cells, including PC, BC, SC, and GC in the forestomachs; for example, ROMO1, which encodes a protein with antimicrobial activity against a variety of bacteria [39] and S100A14, one of the S100A gene families associated with the immune system [40]. Genes related to defense against fungi were highly expressed in the SC and GC; this includes LOC101103771 (also known as S100A12), which exhibits antifungal activity through zinc sequestration [41, 42]. However, the CD40 that defense to protozoan was lowly expressed in all cell types. Genes related to detection, response, and regulation of defense responses to viruses were highly expressed in a wide range of cell types. However, genes related to the responses to microbes, including bacteria, fungi, protozoans, and viruses, were expressed in a cell type-specific fashion in abomasum cells; for instance, S100A14 was highly expressed only in PMC2 and FC. Collectively, we identified cell types with expression of specific genes associated with the host responses to various microorganisms.

Non-immune cells recruit immune cells through MIF signaling

To understand the communication between immune and non-immune cells in coordinating immune responses, we used CellChat to construct ligand-receptor maps. Notably, our analysis revealed that T cells and monocytes/macrophages interacted with non-immune cells through the macrophage migration inhibitory factor (MIF) signaling pathway in the rumen, reticulum, and omasum (Fig. 6A–C). Contribution analysis of ligand-receptor pairs identified that MIF served as the ligand, with CD74/CXCR4 or ACKR3 serving as the receptors (Fig. 6D–F). Notably, the MIF-(CD74/CXCR4) ligand-receptor pair exhibited the highest contribution to MIF signaling in the rumen and omasum (Fig. 6D and F), while MIF-ACKR3 ligand-receptor pair showed the highest contribution in the reticulum (Fig. 6E). Furthermore, we investigated the expression levels of these ligands and receptors in each cell subtype of the forestomach mucosa. Except for granular cell 3 in the rumen and undefined epithelial cells in the omasum, all other cell subtypes expressed high levels of MIF (Fig. 6G–I). Monocytes/macrophages are the major source of CD74, whereas T cells are the major source of CXCR4 and ACKR3. Together, our results revealed that forestomach cells showed a generally high level of MIF, which potentially recruits immune cells into the tissue.

Fig. 6
figure 6

Communication between non-immune cells and immune cells. A–C MIF signaling networks in rumen (A), reticulum (B), and omasum (C), respectively. In the networks, edge width represents the communication probability. D–F The bar plot showing the relative contribution of each ligand-receptor pair to the overall communication MIF signaling networks in the rumen (D), reticulum (E), and omasum (F), respectively. G–I Violin plots showing the expression distribution of ligand and receptor genes of MIF signaling in different cell types, in the rumen (G), reticulum (H), and omasum (I), respectively



Source link