A cell transcriptomic profile provides insights into adipocytes of porcine mammary gland across development

Abstract

Background: Studying the composition and developmental mechanisms in mammary gland is crucial for healthy growth of newborns. The mammary gland is inherently heterogeneous, and its physiological function dependents on the gene expression of multiple cell types. Most studies focused on epithelial cells, disregarding the role of neighboring adipocytes.

Results: Here, we constructed the largest transcriptomic dataset of porcine mammary gland cells thus far. The dataset captured 126,829 high-quality nuclei from physiological mammary glands across five developmental stages (d 90 of gestation, G90; d 0 after lactation, L0; d 20 after lactation, L20; 2 d post natural involution, PI2; 7 d post natural involution, PI7). Seven cell types were identified, including epithelial cells, adipocytes, endothelial cells, fibroblasts cells, immune cells, myoepithelial cells and precursor cells. Our data indicate that mammary glands at different developmental stages have distinct phenotypic and transcriptional signatures. During late gestation (G90), the differentiation and proliferation of adipocytes were inhibited. Meanwhile, partly epithelial cells were completely differentiated. Pseudo-time analysis showed that epithelial cells undergo three stages to achieve lactation, including cellular differentiation, hormone sensing, and metabolic activation. During lactation (L0 and L20), adipocytes area accounts for less than 0.5% of mammary glands. To maintain their own survival, the adipocyte exhibited a poorly differentiated state and a proliferative capacity. Epithelial cells initiate lactation upon hormonal stimulation. After fulfilling lactation mission, their undergo physiological death under high intensity lactation. Interestingly, the physiological dead cells seem to be actively cleared by immune cells via CCL21-ACKR4 pathway. This biological process may be an important mechanism for maintaining homeostasis of the mammary gland. During natural involution (PI2 and PI7), epithelial cell populations dedifferentiate into mesenchymal stem cells to maintain the lactation potential of mammary glands for the next lactation cycle.

Conclusion: The molecular mechanisms of dedifferentiation, proliferation and redifferentiation of adipocytes and epithelial cells were revealed from late pregnancy to natural involution. This cell transcriptomic profile constitutes an essential reference for future studies in the development and remodeling of the mammary gland at different stages.

Keywords: Adipocytes; Cell–cell interaction; Development; Mammary gland; snRNA-seq.

Fig. 1

Generation of a stage-specific single-cell atlas of the pig mammary gland. A Histological observations of the mammary gland at five developmental stages. B Statistical analysis of adipocytes area and size. C Schematic workflow for snRNA-seq sequencing. D UMAP visualization of all clusters colored by all cell types. Seven cell clusters were identified in the dataset. E The DEGs analysis shows upregulated genes (Adjusted P value < 0.01) across all seven clusters. F Cell type annotation for all clusters is provided in the bubble chart. G GO annotation and KEGG pathway analysis of differentially expressed genes in each cell-type. H UMAP illustration of cells colored by clusters in separate development stage. I Stacked bar plots represent the proportions of nuclei in the mammary gland

Fig. 2

Gene expression patterns of epithelial cells at different developmental stages. A Proportional Area Chart (Half Circle). Two groups of half circles indicate two DEGs sets of epithelial cells, and the areas represent the number of DEGs. B Fuzzy clustering of expression data at five developmental stages. Purple or red colored lines correspond to genes with high membership value, and y axis represents the normalized expression value from the Mfuzz result. C Monocle trajectory inference traces a path of pesudotime and group types. D The heatmap reveals the relative gene expression level of 3 clusters at 2 branches based on branched expression analysis modeling, combined with the GO/KEGG enriched items for each cluster

Fig. 3

Subtype classification of epithelial cells in G90 (A) and PI2 (E). Monocle trajectory inference traces a path of pseudotime of epithelial cells in G90 (B) and PI2 (F). The heatmap reveals the relative gene expression of 3 clusters at 2 branches based on branched expression analysis modeling in G90 (C) and PI2 (G), combined with the GO/KEGG enriched items of each cluster. Visualization of the transition of highly expressed genes in pseudotime ordering of epithelial cells in G90 (D) and PI2 (H)

Fig. 4

Gene expression patterns of adipocytes at different developmental stages. A Proportional Area Chart (Half Circle). Two groups of half circles indicate two DEGs sets of adipocytes, and the areas represent the number of DEGs. B Fuzzy clustering of expression data at five developmental points. Purple or red colored lines correspond to genes with high membership value, and y axis represents the normalized expression value from the Mfuzz result. C–D The heatmap reveals the relative gene expression level of adipocytes in G90 (C) and PI2 (D), combined with the GO/KEGG enriched items of each cluster. E Visualization of the transition of highly expressed genes in pseudotime ordering of adipocytes at different developmental stages

Fig. 5

The typical growth factor type ligand-receptor interactions predicted by iTALK between any two cell types. The network plot showed ligand-receptor interactions detected between each two different cell types. In the network, every node showed a cell type, and the thickness of the arrow lines represented the number of ligand-receptor interactions. The arrows labeled the forward (from signaling cell to target cell) and backward signals. The circos plot displayed the names of each ligand-receptor gene pair and the direction. The outside ring of circos plot exhibited cell types, and the inside ring of circos plot exhibited the details of each interaction ligand-receptor pair. The lines inside the circos plot indicated the relative signal strength of the ligand and receptor. The arrow indicated the receptor

Fig. 6

Spatial transcriptome profiles of the mammary gland at G90 and PI2. The epithelial cells and adipocytes were annotated according to gene-makers. Spatial plots showing the expressions of CSN1S1, CSN2, EGFR and ADIPOQ genes

Fig. 7

The immune cellular composition of colostrum and mature milk. A UMAP plot of five clusters from all sequenced milk cells. B Maker genes of the five cell types. C UMAP plot of endothelial cells, epithelial cells, macrophages, monocytes and T cells in two samples. D The proportion of bar plot of five clusters originating from two samples