A cell transcriptomic profile provides insights into adipocytes of porcine mammary gland across development
Journal of Animal Science and Biotechnology volume 14, Article number: 126 (2023)
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.
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.
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.
The mammary gland is an exocrine gland of ectodermal origin, which develops to produce milk for the nourishment of offspring. As a highly dynamic organ, the mammary gland undergoes a limited embryonic development followed by extensive postnatal pubertal development, and further differentiation and tissue remodeling during pregnancy and lactation. The mammary gland is inherently heterogeneous, with various cells performing different functions. The mammary epithelium contains two major cell types, luminal and myoepithelial cells [1, 2], and they are the main components of acini. The lumen of each acinus is hollow with milk secretions during lactation. Milk ejection relies on the contractility of myoepithelial cells . Fibroblasts actively participate in tissue remodeling, synthesizing and secreting collagen, and organizing into bundles in the developing mammary gland . Immune cells participate in and influence branching morphogenesis . Another type of cell, adipocytes, comprise a large portion of the stromal compartment in the adult non-lactating mammary gland, but gradually disappear during pregnancy, freeing up room for the expanding mammary glands . In 2013, Gregor et al.  selectively knocked out X-box binding protein 1 (Xbp1) of adipocyte in mice mammary gland during lactation, causing adipocytes proliferation and lower milk production. In 2014, Vapola et al.  knocked out the peroxisomal membrane protein 2 (Pxmp2) gene in adipocytes of the mice mammary gland, which restricts epithelial cells development and duct formation during pregnancy. In 2018, Wang et al.  reported that adipocytes from mammary gland can become PDGFRα+ preadipocytes and fibroblast-like preadipocytes through dedifferentiation during lactation. In addition, the gene expression level is significantly different between the dedifferentiated cells and the adipocytes in the non-lactating mammary gland. During involution, PDGFRα+ preadipocytes will proliferate and differentiate into adipocytes. Although adipocytes play such a critical role in gland development, most studies focused on its epithelial component, leaving the role of the neighboring adipocytes largely unexplored in both physiologic and pathologic conditions .
Cell fate decisions are largely based on gene transcription. Therefore, it is particularly critical to identity the cell-types and the gene transcription profile of individual cells to understand mammary gland physiology. Recently, single-cell RNA sequencing (scRNA-seq) emerged as a powerful technique to study complex biological systems at single-cell resolution. In 2017, Bach et al.  systematically constructed cell transcriptomic atlas of the epithelial cells from mouse mammary gland across four developmental stages based on scRNA-seq data. The atlas indicated that mammary epithelial cells are not a group that performs a single function, and could be divided into four classes: basal, mature luminal, luminal progenitor, and luminal intermediate cells. Subsequent studies identified immune cells, fibroblasts and endothelial cells in the mammary gland [11,12,13,14]. Individual cells display state-specific expression patterns. For instance, principal component analysis (PCA) exposed gene transcription differences in both basal and luminal cells in pregnant and nonpregnant mice mammary gland . Moreover, Brugge and co-workers confirmed age-dependent alterations in gene expression by analyzing epithelial, stromal, and immune cells in mice mammary gland . Remarkably, mammary gland development and function depends on intricate interactions of the functional epithelial cells with local stromal cells [16, 17]. In other words, clarification of cell–cell interaction is required.
Considering the lack of focus on adipocytes in mammary glands and the limitations of scRNA-seq, it is not surprising they are rarely identified in the mammary gland. In detail, cells over 50 μm in diameter are difficult to capture by microfluidic droplet generators, hindering gene transcription profiling of tissue-derived adipocytes by scRNA-seq approaches. Furthermore, not all studies identify the same mammary cell types and most agree that cell subpopulations such as secretory alveolar cells during lactation have been incompletely profiled due to technical difficulties to isolate intact cells during dissociation . Notably, single-nucleus RNA-seq (snRNA-seq) has become instrumental to interrogate oversized cells or in complex tissues that are not easily dissociated. This approach enabled the comprehensive mapping of mammary cell transcription in different physiological states. Regrettably, only one report identified adipocytes in mammary glands of adult human using snRNA-seq as of June 2023 , which means that, to date, no gene transcription data from adipocytes of porcine mammary glands is available.
The physiological function of a tissue is not only dependent on the gene expression of its individual cells. The internal spatial organization of these cells is also critically important . However, spatial information of individual cells is lost during tissue dissociation for snRNA-seq sequencing [20, 21]. Contrastingly, spatial transcriptomics (ST), a more recent method, enables the visualization and quantitation of the transcriptome in individual tissue sections, retaining spatial molecular information . However, the spatial expression pattern of adipocytes in healthy porcine mammary glands has not been reported. Thus, we describe the spatial transcriptome profiles of adipocytes in non-lactation mammary gland. Furthermore, this spatial expression profile validates the accuracy of the cell type annotation based on snRNA-seq data.
In this study, we constructed the largest transcriptomic dataset of the porcine mammary gland thus far to reveal its development at 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 in this dataset, including epithelial cells, adipocytes, endothelial cells, fibroblasts cells, immune cells, myoepithelial cells and precursor cells. Their gene transcription patterns were determined, and an interaction network between adipocytes and other cell groups was constructed at both snRNA-seq and spatial transcriptomics levels. To explore the impact of developmental stages on milk composition, we identified the type and proportion of milk secreted during colostrum and mature milk using scRNA-seq. Cell phenotype analysis showed more abundant macrophages in colostrum than mature milk, which likely explains why colostrum possess innate immune activity. The swine cell atlas here reported will guide future studies in mammary physiology. It also provides insights into the development stage-specific gene expression profiles of epithelial cells and adipocytes.
Materials and methods
Sample collection and histological observation
Multiparous sows in second breeding cycle were divided into 5 groups [d 90 of gestation (G90); d 0 (L0), d 20 (L20) after lactation; 2 d (PI2) and 7 d (PI7) post natural involution]. Each group has 4 independent biological replicates. First parity sows usually face higher risk and stress response. Thus, we selected the multiparous sows as an experimental animal to ensure the reliability and reproducibility of the results. All pigs were fed well-characterized normal diets, according to the nutritional requirements outlined by the Feeding Standard of Swine (NY/T 65–2004) and published by the Ministry of Agriculture and Rural Affairs of the People’s Republic of China . Sows were humanely euthanized at five developmental stages for collection of the mammary glands (approximately 3 cm3) from the third region on the right side (Fig. 1A). Part of fresh tissue was used to prepare frozen sections for histological observation (HE staining, n = 4) (G90, L0, L20, PI2 and PI7) and spatial transcriptome sequencing (n = 1) (G90 and PI2). Meanwhile, the remaining tissue (n = 2) (G90, L0, L20, PI2 and PI7) was collected, immediately snap-frozen in liquid nitrogen, and then transferred to −80 °C until further nuclei isolation. In addition, 15 mL fresh colostrum (L0) and mature milk (L20) were collected for scRNA-seq (n = 1) for cellular component analysis.
Single nuclei/cell RNA sequencing
The snRNA-seq and scRNA-seq were respectively performed for the construction of the transcriptional atlas of the mammary glands and milk cellular components. Nuclei isolation was carried out using GEXSCOPE® Nucleus Separation Solution (Singleron Biotechnologies, Nanjing, China) according to the manufacturer's product manual. Isolated nuclei were resuspended in PBSE to 106 nuclei per 400 μL, filtered through a 40-μm cell strainer, and counted with Trypan blue. The concentration of single nuclei suspension was adjusted to 3–4 × 105 nuclei/mL in PBS. Subsequently, single nuclei/cell suspension was loaded onto a microfluidic chip (GEXSCOPE® Single NucleusRNA-seq Kit/GEXSCOPE™ Single-Cell RNA Library Kit, Singleron Biotechnologies) and snRNA-seq/scRNA-seq libraries were constructed according to the manufacturer’s instructions (Singleron Biotechnologies). Finally, the resulting snRNA-seq/scRNA-seq libraries were sequenced on an Illumina HiSeq × 10 instrument with 150 bp paired end reads.
Data processing and cell-type annotation
Sequencing reads were processed using the CeleScope v1.1.7 pipeline (Singleron). Then, raw reads were aligned to the Sscrofa11.1 reference genome, generating count matrices. To remove low-quality droplets, we excluded any nuclei or cells expressing less than 300 genes and more than 25% mitochondrial genes. After filtering, 126,829 nuclei were reserved and used for dimension-reduction and clustering. The normalization and scale of all gene expression values were performed by NormalizeData and ScaleData function. Principal component analysis (PCA) relied on the top 2000 variable genes screed by FindVariableFeatures. Next, integration of snRNA-seq and scRNA-seq data was respectively done using the harmony package (https://github.com/immunogenomics/harmony) to control for batch effects when integrating data from different development samples. For dimensionality reduction, the top 15 principal components were selected to calculate 2D dimensional reductions by Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) for all sequencing libraries based on an elbow approach. Finally, cell clusters were identified using the Louvain algorithm at a resolution of 0.3, implemented by the FindCluster function of Seurat (v4.3.0).
The clusters were partitioned into distinct cell types and annotated by the expression of known marker genes. The expression signatures of cell‐type‐specific genes were detected using the “FindAllMarkers” function. The criteria to identify cell‐type‐specific genes were set as follows: absolute log2 fold change (FC) ≥ 1 and the minimum cell population fraction in either of the two populations was 0.25. The expression pattern of each marker gene was visualized by applying the “DotPlot” function in Seurat.
Gene function analysis
To reveal gene transcription dynamics of adipocytes and epithelial cells, we detect time-series expression profiles of differentially expressed genes (DEGs) between adjacent development stages using Mfuzz algorithm. DEGs were identified using FindMarkers function in Seurat (log2FC > 1, P < 0.05). The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichments were conducted to identify the main function of the DEGs having the same expression trend using the ClusterProfiler v4.0.2 package.
In addition, we identified DEGs (log2FC > 1, P < 0.05) between any two adjacent developmental stages of the annotated seven cell types. To clarify their biological function, GO annotation and KEGG pathway analysis were performed, generating 28 sets of significant GO terms and KEGG pathways. We compared four sets of significant terms and pathways in each cell type, and visualized the results using Venn Diagrams. Another group of Venn diagrams displayed the intersection of GO terms or KEGG pathways and were annotated in the seven type cells in different physiological processes.
To characterize the physiological state of epithelial cells and adipocytes in G90 and PI2, we calculate their differentiation trajectory based on default parameters of DDRTree method in Monocle2 package. In detail, normalized unique molecular identifier (UMI) count was fed as the input for Monocle2. Genes with a residual greater than 1× the estimated mean–variance split, were identified as high-dispersion genes using the “estimateDispersions” method. After running the “setOrderingFilter” function, dimensionality reduction was applied to the data with the default parameter of DDRTree method . The trajectory was visualized by plot_cell_trajectory function. Furthermore, branch analysis was performed by branched expression analysis modeling (BEAM), and visualized via the plot_genes_branched_heatmap function .
The high-confidence ligand (L)—receptor (R) interactions were performed to investigate the interaction between mammary gland cells in five development stages using the iTALK package . Expressed genes were selected for L-R interaction analysis according to the following criteria: (1) the top 20 highly expressed genes and (2) marker genes in corresponding clusters.
Spatial transcriptome sequencing of mammary gland
Spatial expression of mammary gland during G90 and PI2 was performed to evaluate cell-type annotation accuracy. In detail, the 10 μm frozen tissue sections were placed on one of the Visium gene expression slide capture areas in a slide. The RNA quality of mammary glands was assessed by Agilent 2100, and RNA integrity number (RIN) of tissues greater than 7 were used for Visium spatial gene expression experiments. The Visium Spatial Gene Expression Slide & Reagent kit (10 × Genomics) was used to construct sequencing libraries according to the Visium Spatial Gene Expression User Guide (CG000239, 10 × Genomics). Tissue permeabilization was optimized during the tissue optimization procedure. Reverse transcription experiment and sequencing libraries were then prepared following the manufacturer’s protocol. Sequencing was performed with a Novaseq PE150 platform according to the manufacturer’s instructions (Illumina) at an average depth of 300 million read-pairs per sample.
We used in-house script to perform basic statistics of raw data, and evaluate the data quality and GC content along the sequencing cycles. Raw FASTQ files and histology images were processed by sample with the Space Ranger (version spaceranger-1.2.0, 10 × Genomics) software with default parameters. The filtered gene-spots matrix and the fiducial-aligned low-resolution image was used for down-streaming data analyses.
Reads were demultiplexed using Space Ranger software v.1.0.0 (10 × Genomics) and annotated with the reference genome Sscrofa11.1. Subsequently, the generated count matrices were loaded into Seurat environment, and the data were normalized (using the ‘SCTransform’ function in Seurat for independent tissue sections), reduced and visualized. SPOTlight analysis was performed for deconvolution analysis as Elosua-Bayes et al. reported . In brief, the proportion of signature of the selected snRNA-seq cell type is equal to the sum of the proportions of each cell type in different regions, divided by the sum of the proportions of that cell type in all spots.
The statistical analysis of the area and diameter of adipocytes were performed using GraphPad Prism 9 (GraphPad, San Diego, CA, USA) and one-way ANOVA. Data were presented as the mean ± standard deviation (SD), and P < 0.05 (*) indicated a significant difference.
Sequencing data were deposited in the Gene Expression Omnibus (GEO) with the accession code GSE227425.
Generation of a development stage-specific single-cell atlas of mammary gland
We collected mammary glands of ten female pigs across five developmental stages, including d 90 of gestation (G90); d 0 (L0) and 20 (L20) after lactation; 2 d (PI2) and 7 d (PI7) post natural involution for histological observations (Fig. 1A). Mammary gland sections stained with HE indicated that acini appeared in late gestation (G90). In addition, the mammary gland went through remarkable morpho-functional changes in its adipocytic components, i.e., both area and size of adipocytes were greater in non-lactation than lactation stages (Fig. 1B and Table S1).
To understand the tissue composition and gene-expression dynamics, we generated snRNA-seq profiles from 10 mammary glands across the five developmental stages (Fig. 1C). After quality control, a transcriptomic dataset with 126,829 high-quality nuclei were retained, with numbers ranging from 18,693 in the d 90 gestation group to 9,331 in the 2 d post natural involution group (Fig. S1 and Table S2). On average, we detected 2,542 UMIs and 1,395 genes per nuclei. Upon batch effect correction, the 126,829 nuclei separated into multiple clusters (Fig. 1D) using UMAP . Seven cell clusters were defined according to the expression levels of specific markers (Fig. 1E–F and Table S3) including adipocytes, epithelial, fibroblasts, endothelial, myoepithelial, immune and precursor cells (Fig. 1D). For instance, adiponectin (ADIPOQ) and epithelial cell adhesion molecule (EPCAM) expression subdivide adipocytes and epithelial cells from mammary glands, respectively. The assigned cell types were further confirmed by function analysis of gene sets identified by “FindAllMakers” function in Seurat  (Fig. 1G). Unsurprisingly, cellular functions derived from shared gene annotations were associated with phenotypic similarity. The biological functions of the gene sets were linked with the cell types, such as the gland development term, the epithelial cell development term and the lactation term which were only significantly enriched in the annotated epithelial cell. Similarly, the regulation of lipolysis in adipocytes pathway, the PPAR signaling pathway and the adipocytokine signaling pathway were significantly enriched in the annotated adipocytes. This provided further evidence of cell type identification accuracy.
The total number of each cell type ranged from 61,211 (48.26%) for epithelial cells, to 687 (0.54%) for adipocytes in the mammary gland (Fig. 1H and Table S4). In addition, cell-type composition dynamics changed in mammary glands during different developmental stages (Fig. 1I). Globally, adipocytes constituted between 0.34% and 0.09% of all cells present, with a higher proportion in non-lactation period compared with lactation. This is consistent with our phenotypic results (Fig. 1B) and previous study . Conversely, the proportion of epithelial cells, varying from 72.43% to 19.37%, was higher in lactation versus non-lactation period, gradually decreasing gradually with mammary gland remodeling. In addition, the ratio of fibroblasts, endothelial and myoepithelial cells was greater in non-lactation period rather than lactation.
Gene expression patterns of epithelial cells across five developmental stages
To unveil the gene transcription dynamics of epithelial cells across adjacent developmental stages, we identified differentially expressed genes (DEGs) using the FindMarkers function in Seurat (Fig. 2A and Table S5). The time-series gene expression profiling of epithelial cells exhibits six time-dependent expression patterns (Fig. 2B). During late pregnancy (G90), the main functions of highly expressed genes (cluster 5) are regulation of cell division and cell differentiation. After deliver (L0), specific gene-expression (cluster 6) regulated by hormone stimulation initiate lactation. As the newborn grows (L20), the biosynthetic process (cluster 2) is activated in maternal mammary epithelial cells. Upon cessation of suckling by the offspring, the involution of the mammary gland is initiated. In early natural involution (PI2), most of highly expressed genes in epithelial cells (cluster 3) participated in apoptotic and cell proliferation inhibition to allow space for other cell types. At 7 d post natural involution (PI7), the immune mechanism is activated.
Although epithelial cells are captured from the same developmental stage, they are still in different physiological states. Monocle provides a convenient way to screen all pseudotime-dependent genes and identify genes following similar kinetic trends. In late pregnancy (G90), the pseudotime analysis shows three cell subsets with similar gene expression patterns (Fig. 2C). After function analysis, cells from the three subsets performed different functions according to pseudotime. In cell cluster 1, highly expressed genes are mainly enriched in terms or pathways related to cell differentiation, development and morphogenesis. In cell cluster 2, genes participated in hormone signal and the development process of mammary gland, alveolus and lobules. In cell cluster 3, high expression genes promoted lactation and substance synthesis in epithelial cells. Similarly, three cell clusters were identified in early natural involution (PI2) (Fig. 2D). The pseudotime series analysis reveals three physiological states of these cells: cessation of proliferation in response to hormone stimulation (Rap1 signaling pathway), programmed cell death (IL-17 signaling pathway and NF-kappa B signaling pathway), and dedifferentiation into mesenchymal-like cells (mesenchymal cell differentiation).
To characterize the physiological state of epithelial cells in G90 and PI2, we used Monocle2 package to calculate the differentiation trajectory of these cells. According to known marker genes [31,32,33], cells were classified as epithelial precursor cells (ALDH1A3, CD14 and KIT), luminal epithelial cells (CSN2 and LALBA) and hormone-sensing epithelial cells (ESR1, PRLR and PGR) (Fig. 3A and Fig. 3E). Cell differentiation trajectories indicate that hormone-sensing cells and luminal cells originate from epithelial precursor cells (Fig. 3B and Fig. 3F).
In late pregnancy (G90), gene function analysis revealed hormone receptor genes are activated in hormone-sensing epithelial cells, including erb-b2 receptor tyrosine kinase 4 (ERBB4), transforming growth factor alpha (TGFA) and estrogen receptor 1 (ESR1) (Fig. 3C–D). Functional activation of luminal cells depended on Acetyl-CoA carboxylase (ACACA), fatty acid synthase (FASN) and ACLY (Fig. 3C–D).
In early natural involution (PI2), the hormone-sensing epithelial cells received hormone stimulation [phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1), mitogen-activated protein kinase 10 (MAPK10) and TGFA] and activate apoptosis related terms and pathways (BMPR1A, BMPR1B and STK3). At the same time, luminal epithelial cells activated the positive regulation of programmed cell death term [superoxide dismutase 2 (SOD2), cathepsin C (CTSC) and beta-1,4-galactosyltransferase 1 (B4GALT1) and the regulation of epithelial cell apoptotic process pathway (beta-2-microglobulin (B2M), nuclear protein 1, transcriptional regulator (NUPR1) and programmed cell death 4 (PDCD4)] (Fig. 3G–H).
Gene expression patterns of adipocytes across five developmental stages
To unveil gene transcription dynamics of adipocytes across adjacent developmental stages, we identify DEGs using the FindMarkers function in Seurat (Fig. 4A and Table S5). Based on DEGs identified from adipocytes, we employed Mfuzz algorithm to detect time-series expression profiles of genes across developmental stages. Six time-dependent expression patterns were characterized in adipocytes and investigated for their biological significance (Fig. 4B). The DEGs of adipocytes in cluster 1 and 5 show high expression levels during non-lactation. The function of these genes is significantly enriched in cellular development and macromolecule metabolic processes term. This is consistent with the results of both our phenotypic data (Fig. 1A–B) and previous research . The DEGs in cluster 2 are specific expression at colostrum stage (L0) versus the other four stages. These DEGs are mainly involved in the response to hormone and adipocytokine. In cluster 3 and cluster 4, the DEGs are highly expressed at L20, and sustain adipocytes survival. This correlates a compression of the adipocyte space during lactation. Highly expressed genes in cluster 6 are mostly involved in the response to hormone stimulus and cell proliferation during early natural involution (PI2).
In many biological processes, cell growth, differentiation and development do not progress in perfect synchrony. Single-cell expression studies of cell differentiation often capture cells distributed across the entire process. To understand the gene gradient expression at different developmental stages, we applied pseudotime and trajectory analysis for adipocytes using Monocle 2. In late pregnancy (G90), the majority of adipocytes are in a highly differentiated state. When the living space is occupied by other cell types (Fig. 1A), adipocytes gradually dedifferentiate during the onset and maintenance of lactation (Fig. 4C). During lactation, highly expressed genes regulate cell morphogenesis and keep adipocytes at undifferentiated state. These highly expressed genes are significantly enriched in the canonical insulin receptor signaling pathway and AMPK signaling pathway (Fig. 4D). In the two pathways, we found three key genes, ZFP36 ring finger protein like 1 (ZFP36L1), forkhead box O1 (FOXO1), and lipid phosphate phosphohydrolase 1 (LPIN1).
In natural involution, the highly expressed genes regulate adipocytes redifferentiation and lipid metabolism (the fat cell differentiation term and the triglyceride metabolic process term) (Fig. 4D). The pseudo time series analysis indicated that adipocytes displayed a mature state at the PI2 and PI7 stages, with terms related to adipocyte substance synthesis, including the regulation of lipid storage term (CD36 molecule, CD36), the response to fatty acid term (lipoprotein lipase, LPL), and the triglyceride metabolic process term (ATP binding cassette subfamily A member 1, ABCA1) (Fig. 4E).
Many cell types of the mammary gland contribute to its structure, development, and ultimate function in a dynamic and reciprocal fashion. As proof of principle of the application of this dataset for describing mammary gland cell–cell interactions, we next detected the distribution of cytokine, growth factors and other receptors across the seven type cells. The number of paired ligand-receptor (L-R) interactions showed in the network plots reveals the interactions between each two different cell types and within the same cell type. The Circos plots display the top 20 L-R pairs (Fig. 5 and Fig. S2). In terms of growth factors, the interactions between epithelial cells and other cells are most common during d 90 gestation than the other four stages. Aside from adipocytes, epithelial cells strongly bind endogenous and exogenous growth factor ligands secreted by other cell types. Interactive networks displayed the L-R pairs of growth factors which promoted epithelial cell development, these were mainly driven by TGFB2-TGFBR3, TGFA-EGFR and TGFBR1-TGFBR3 pathways (Fig. 5 and Fig. S2). At the same time, endothelial cells promote the growth of other five cell types, except adipocytes, by secreting TGFB2 and PDGFD ligand. Meanwhile adipocytes showed weak interactions with other cells. Interestingly, the top 20 cell interactions suggest that the interaction between precursor cells and other cell populations only occurs in the transition state of non-lactation and lactation, that is, G90 and PI2 periods (Fig. 5A and D).
At the onset of lactation (L0), the cytokine-type L-R pairs (CXCL12-ITGB1) connects endothelial cells with other cells, except adipocytes and precursor cells. Another important feature of the L0 period, is that the epithelial, myoepithelial and immune cells communicate with endothelial cells via PDGFC-KDR pathway (Fig. 5 and Fig. S2). Fibroblasts secret VEGFC ligand to combined with FLT1, FLT4, LYVE1, ITGA9 and ITGB1 receptors from endothelial cells. At mature milk stage (L20), the cell–cell interactions show that the cytokine ligands secreted by endothelial cells (CCL5, CCL21) and myoepithelial cells (CCL2) activate the adipocyte surface receptor ACKR4 (Fig. 5C and Fig. S2). In addition, two autocrine pathways (FGF10-FGFR1 and ADIPOQ-ADIPOR2) were activated during this period.
In 2 d post natural involution (PI2), adipocytes frequently received cytokine signals from other cell types (Fig. 5D and Fig. S2). Notably, different from mature milk stage, the autocrine pathways maintained the adipocyte proliferation through FGF1-EGFR and IGFBP4-LRP6 pairs at PI2, rather than FGF10-FGFR1 pairs (Fig. 5 and Fig. S2). In epithelial cells, the positive regulation of programmed cell death term was activated. All mammary cells sent out the demand for proliferation and differentiation to immune cells via IL34-CSF1R pairs. In post natural involution (PI2 and PI7), epithelial cells stimulated myoepithelial cell growth through the growth factor type ligand-receptor (TGFA-EGFR) pathway during early natural involution (Fig. 5 and Fig. S2). In PI7, TGFA regulated the proliferation of epithelial cells through an autocrine signaling pathway (TGFA-ERBB4).
Epithelial cells and adipocytes annotated in situ with precise spatial resolution
Mature acini, the basic unit of galactosis, are observed in the G90, whose main component is mammary epithelial cells. Acinar maturation indicates that the mammary gland has been fully prepared for lactation during G90. Additionally, casein alpha s1 (CSN1S1) and casein beta (CSN2) are the key genes for mammary epithelial cells participation in lactation. However, snRNA-seq loses spatial information of the expression profiles of CSN1S1 and CSN2. Thus, we verified that the two genes were specifically expressed in the acinar region using spatial transcription sequencing (Fig. 6A and Fig. S3). In detail, we integrated the snRNA-seq and spatial datasets based on SPOTlight with a deconvolutional procedure, generating the spatial distribution of sequenced cells. The spots with CSN1S1 and CSN2-specific expression are largely overlapping with the region annotated as epithelial cells. We found that EGFR expression was a specific expression marker in epithelial cells from G90 (Fig. 6A), and repopulation of adipocytes a distinctive feature of natural involution (PI2) (Fig. 1A). The spatial distribution of adipocytes annotated by SPOTlight showed a consistent spatial colocalization of adipocytes in HE-stained section. ADIPOQ is a marker of adipocytes that encodes a protein hormone, adiponectin, involved in the regulation and inhibition of lipogenesis and the stimulation of fatty acid oxidation. Its spatial expression pattern also supports the annotation of snRNA-seq data (Fig. 6B and Fig. S3).
Single-cell sequencing of milk cells from colostrum and mature milk
The lactation stage affects the cellular component of milk . To explore this effect, we collected fresh colostrum (L0) and mature milk (L20) for scRNA-seq. A total of 31,943 high-quality cells were captured from colostrum and mature milk, including endothelial cells, epithelial cells, macrophages, monocytes and T cells (Fig. 7A). The cell-type classification showed a higher number of macrophages in colostrum (12.06%) versus mature milk (9.70%) (Fig. 7C and D).
The mammary gland provides essential nutrients to the suckling infant [35, 36]. In most mammals, mammary gland morphogenesis begins at embryonic period. After birth, it undergoes three successive stages: puberty, pregnancy and lactation, and natural involution [31, 37, 38]. To date, the adipocyte’s transcriptome profile in swine mammary glands has not been reported. Therefore, the mechanisms by which adipocytes regulate porcine mammary developments and participate in remodeling remains to be elucidated. Here, we measured the gene expression of swine mammary glands across five development stages using snRNA-seq, capturing porcine mammary adipocytes. In this dataset, another six cell types were also annotated; epithelial, fibroblasts, endothelial, myoepithelial, immune and precursor cells. In addition, histological observations revealed a dynamic change in the area and size of the adipocytes through mammary gland development, and indicated that acini appeared in late gestation (G90), which implies cells in mammary gland are preparing for lactation at this developmental stage [39, 40].
In late pregnancy (G90), the majority of adipocytes are highly differentiated, with minimal communications with other cells. Such minimal communications means that adipocytes are less regulated by growth factors than the other six cell groups in G90. We believe this is the reason for the lower adipocyte proportion in the subsequent L0 period. Functional analysis of gene clusters that control cell fate unveiled the activation of GO terms negatively regulated to cell morphogenesis and cell proliferation. This may be the driving force for decreased fat cell volume and number during lactation. In the activated terms and pathways, we found a key gene, IGF1R, which regulates cell size through c-Myc family members and plays an important role in the process of cell dedifferentiation, consistent with previous studies [41,42,43].
Another characteristic event of G90 is that the lactation function of epithelial cells is improved, as specifically shown by the activation of genes related to cell division and differentiation. Most epithelial cells undergo three stages to achieve lactation functions, namely, cellular differentiation, hormone sensing, and metabolic activation. Functional gene analysis revealed that hormone-sensing epithelial cells promote mammary gland development by activated hormone receptor genes, such as ERBB4, TGFA and ESR1. ERBB4 is required for the differentiation of mouse mammary epithelial cells during pregnancy and promotes differentiation of murine and human mammary epithelial cells in cell culture . However, other animal models demonstrated that TGFA promotes epithelium growth in mammary glands by binding to the EGF receptor, activating its kinase cell signaling. Meanwhile, luminal cells perform essential functions in milk synthesis and secretion (Fig. 3G and H). ACACA, FASN and ACLY are key genes to improving the lactation function of luminal epithelial cells. ACACA and FASN are two critical genes required for fatty acid synthesis in milk, mainly acting in the elongation of the fatty acid chain [45, 46]. ACLY converts cytoplasmic citrate to acetyl-CoA and oxaloacetate and catalyzes the first step of the de novo lipogenesis pathway . Epithelial cell maturity is marked by higher expression levels of genes associated with casein synthesis [48, 49], such as CSN1S1 and CSN2, whose expression are detected at single cell and spatial levels.
At the onset of lactation (L0), the proliferation of mammary cells, except adipocytes and precursor cells, is strongly affected by the cytokine-type L-R pairs, CXCL12-ITGB1. The team of Ryota Kawahara proved that the phenotype of ITGB1-KO mice is embryonic lethal . In L0, the epithelial, myoepithelial and immune cells all promote growth, proliferation, migration of endothelial cells and vascular network formation through PDGFC-KDR pathway (Fig. 5 and Fig. S2) . Fibroblasts regulate endothelial cell differentiation and promote angiogenesis by secreting VEGFC ligand, one of the strongest modulators of angiogenesis , via FLT1, FLT4, LYVE1, ITGA9 and ITGB1 receptors [53,54,55,56,57]. These L-R pathways jointly promote angiogenesis and ensure the transportation of nutrients required for lactation [58, 59].
Epithelial cells initiate lactation (L0) through hormone signals. As lactation progresses, epithelial cells actively synthesize milk to meet nutritional requirements for newborns. After fulfilling lactation mission, epithelial cells undergo programed cell death (L20). Cell–cell interaction analysis indicated that epithelial cells promoted the proliferation and differentiation of immune cells by secreting IL34 ligand and binding to CSF1R surface receptor [60, 61]. This suggests that epithelial cells can be actively cleared by immune cells to eliminate dead cells under high intensity lactation. This biological process may be an important mechanism for maintaining homeostasis of the mammary gland. It has been established that immune cells including neutrophils, macrophages and lymphocytes are present in milk . Numerous data on animal studies have shown that maternal immune cells can be transferred to newborns through milk. Langel et al.  suggests that these maternal immune cells contribute to the maturation of the innate immune system in the offspring. There are abundant macrophages in colostrum. With the lactation process, the number of macrophages in mature milk decreases [64, 65]. Importantly, these macrophages often elicit significant and non-homeostatic inflammatory responses. In 2020, Zimmermann and Macpherson  demonstrated that T-regulatory cells (Treg) can be transmitted to the offspring via milk, and persist for a long time in a mechanism that provides microbial and pathogenic resistance to the offspring. These results suggest that Treg in milk can help newborns build an early immune barrier.
Our data indicates that the adipocytes area accounts for less than 0.5% during lactation (L0 and L20). This results from the dedifferentiation mechanism of adipocytes. To assure survival in the extremely stressful living environment, the adipocyte maintains a poorly differentiated state with high proliferative capacity, which may be driven by hormone inhibition of adipocyte differentiation (the canonical insulin receptor signaling pathway and AMPK signaling pathway). In the two pathways, three key genes (ZFP36L1, FOXO1 and LPIN1) inhibit the differentiation of adipocytes. ZFP36L1 inhibits intracellular fat synthesis , while FOXO1 gene inhibits adipocyte differentiation . Previous studies have shown that LPIN1 expression is higher in preadipocytes than in mature adipocytes [69, 70]. Similarly, LPIN1 displays a higher expression in L0 and L20 rather than G90. At mature milk stage (L20), the cell–cell interactions show that the cytokine ligands secreted by endothelial cells (CCL5, CCL21) and myoepithelial cells (CCL2) activate the adipocyte surface receptor ACKR4 (Fig. 5C and Fig. S2). After receiving CCL21 signals from other cells, adipocytes reportedly activate ACKR4 to remove chemokines and evade phagocytosis by inflammatory cells [71,72,73]. A previous study confirms that the absence of FGFR1 gene leads to a delay in mammary gland development, with a short-term decrease in cell proliferation . This may indicate that the autocrine FGF10-FGFR1 pathway is the required for proliferation potential of adipocytes during L20 (Fig. 5 and Fig. S2). Besides, adipocytes inhibit apoptosis via the ADIPOQ-ADIPOR2 pathway [6, 75]. Overall, adipocytes are actively mobilizing their growth potential in the period of physiological lactation degradation and preparing for remodeling.
Previous research has confirmed that the mammary gland undergoes the epithelial programmed cell death during lactational involution . Our latest work anatomizes this process in more detail. Epithelial cell populations first stop proliferating (Rap1 signaling pathway), then undergo programmed cell death (NF-kappa B signaling pathway and IL-17 signaling pathway), and the last surviving cells dedifferentiate into mesenchymal stem cells (BMPR1A, TGFBR3 and WWTR1). In detail, the Rap1 signaling pathway plays an important role in cell proliferation , while the NF-kappa B signaling pathway and IL-17 signaling pathway drive apoptosis in epithelial cells [78, 79]. Non-apoptotic epithelial cells are dedifferentiated to mesenchymal cells as “seeds” for subsequent lactation cycles. The transforming growth factor beta receptor 3 (TGFBR3) and WW domain containing transcription regulator 1 (WWTR1) induce the transformation of epithelial cells into mesenchymal cells [80, 81]. Additionally, the bone morphogenetic protein receptor type 1A (BMPR1A) plays an important role in maintaining the undifferentiated state of mammary epithelial cells . Volume reduction is a critical feature in mammary gland involution after weaning. To maintain gland homeostasis, immune cells are required to clear the programmed dead cells . Surprisingly, these cells seem to be actively cleared by immune cells. This conclusion is based on the results of the gene transcription dynamics and cell–cell interaction analysis. At this stage, mammary cells signal the demand for proliferation and differentiation to immune cells via IL34-CSF1R pairs [31, 82]. Similarly, epithelial cells stimulate myoepithelial cell growth through the growth factor type ligand-receptor (TGFA-EGFR) pathway during early natural involution (Fig. 5 and Fig. S2) . We speculate that this dynamic transcription profile is related to epithelial cells role in ejection of milk in the acini through the myoepithelial cell contraction to reduce the inflammatory reaction in the mammary gland. In 7 d post natural involution (PI7), TGFA ligand activates the receptor ERBB4 on the surface of adjacent epithelial cells to maintain proliferation, so that the epithelial cells persist in an extremely narrow living space . This molecular mechanism maintains the lactation potential of the mammary gland and guarantees the next lactation cycle.
Adipocytes area increase by 22.85 percentage points at PI7 versus L20. The regulation of lipid storage and lipid biosynthetic process leads to extracellular lipid transport into the cytoplasm and fusion with other lipid droplets to increase the size of fat adipocytes, which is one of the important ways to hypertrophy . This is consistent with a report from Zwick et al.  that the regulatory mechanism of adipocyte hypertrophy and reoccupation during mammary gland remodeling in mice. In remodeling of mammary glands, the regulation of lipid storage term (CD36), the response to fatty acid term (LPL), and the triglyceride metabolic process term (ABCA1) are activated. The ABCA1 can regulate adipocyte lipid metabolism by adjusting the lipid content of adipose tissue, glucose tolerance, and insulin sensitivity [85, 86]. CD36 is a scavenger receptor that plays a role in adipose energy storage . The LPL gene mainly participates in the uptake of lipids by adipocytes, and is a critical regulatory factor for lipid accumulation in adipocytes, as well as a marker for adipocyte differentiation [88, 89]. In summary, the lineage trajectory analysis revealed that adipocytes underwent dedifferentiation, proliferation and redifferentiation from late pregnancy to natural involution (PI7) [6, 9]. Adipocytes frequently received cytokine signals from other cell types (Fig. 5D and Fig. S2). This phenotype was closely related to the active genes involved in proliferation and differentiation. In PI2, the surface cytokine receptor, ACKR4, participated in adipocyte chemokine’s clearance [77, 78, 90]. Notably, different from mature milk stage, the autocrine pathways maintain adipocyte proliferation through FGF1-EGFR and IGFBP4-LRP6 pairs at PI2, rather than FGF10-FGFR1 pairs (Fig. 5 and Fig. S2) [79, 80]. The above molecular mechanisms maintain the lactation potential of the mammary gland and guarantee the next lactation cycle.
Taken together, this dataset is the largest cell transcriptomic profile of porcine mammary gland across development to date. Herein, we reveal the internal factors of proliferation, differentiation and apoptosis of epithelial cells and adipocytes across five development stages. A vital finding was that epithelial cells are converted into mesenchymal stem cells during the remodeling process. Of note, our data annotated adipocytes in the porcine mammary glands, and clarified the molecular mechanism of dedifferentiation, proliferation and redifferentiation in adipocytes, from late pregnancy to natural involution. Overall, our data provide novel and fundamental insights into the mammary gland development.
Differentially expressed genes
Gene Expression Omnibus
Hematoxylin and Eosin
Kyoto Encyclopedia of Genes and Genomes
Principal component analysis
Single-cell RNA sequencing
Uniform manifold approximation and projection
Al-Khalaf HH, Ghebeh H, Wakil SM, Al-Mohanna F, Aboussekhra A. Interleukin-8 dedifferentiates primary human luminal cells to multipotent stem cells. Mol Cell Biol. 2020;40(9):e00508–19. https://doi.org/10.1128/mcb.00508-19.
Ewald AJ, Brenot A, Duong M, Chan BS, Werb Z. Collective epithelial migration and cell rearrangements drive mammary branching morphogenesis. Dev Cell. 2008;14(4):570–81. https://doi.org/10.1016/j.devcel.2008.03.003.
Duong MT, Akli S, Wei C, Wingate HF, Liu W, Lu Y, et al. LMW-E/CDK2 deregulates acinar morphogenesis, induces tumorigenesis, and associates with the activated b-Raf-ERK1/2-mTOR pathway in breast cancer patients. PLoS Genet. 2012;8(3):e1002538. https://doi.org/10.1371/journal.pgen.1002538.
Li X, Liu R, Su X, Pan Y, Han X, Shao C, et al. Harnessing tumor-associated macrophages as aids for cancer immunotherapy. Mol Cancer. 2019;18(1):177. https://doi.org/10.1186/s12943-019-1102-3.
Aurora AB, Olson EN. Immune modulation of stem cells and regeneration. Cell Stem Cell. 2014;15(1):14–25. https://doi.org/10.1016/j.stem.2014.06.009.
Wang QA, Song A, Chen W, Schwalie PC, Zhang F, Vishvanath L, et al. Reversible de-differentiation of mature white adipocytes into preadipocyte-like precursors during lactation. Cell Metab. 2018;28(2):282–8.E3. https://doi.org/10.1016/j.cmet.2018.05.022.
Gregor MF, Misch ES, Yang L, Hummasti S, Inouye KE, Lee AH, et al. The role of adipocyte XBP1 in metabolic regulation during lactation. Cell Rep. 2013;3(5):1430–9. https://doi.org/10.1016/j.celrep.2013.03.042.
Vapola MH, Rokka A, Sormunen RT, Alhonen L, Schmitz W, Conzelmann E, et al. Peroxisomal membrane channel Pxmp2 in the mammary fat pad is essential for stromal lipid homeostasis and for development of mammary gland epithelium in mice. Dev Biol. 2014;391(1):66–80. https://doi.org/10.1016/j.ydbio.2014.03.022.
Colleluori G, Perugini J, Barbatelli G, Cinti S. Mammary gland adipocytes in lactation cycle, obesity and breast cancer. Rev Endocr Metab Disord. 2021;22(2):241–55. https://doi.org/10.1007/s11154-021-09633-5.
Bach K, Pensa S, Grzelak M, Hadfield J, Adams DJ, Marioni JC, et al. Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA sequencing. Nat Commun. 2017;8:2128. https://doi.org/10.1038/s41467-017-02001-5.
Pal B, Chen Y, Vaillant F, Jamieson P, Gordon L, Rios AC, et al. Construction of developmental lineage relationships in the mouse mammary gland by single-cell RNA profiling. Nat Commun. 2017;8:1627. https://doi.org/10.1038/s41467-017-01560-x.
Li CM, Shapiro H, Tsiobikas C, Selfors LM, Chen H, Rosenbluth J, et al. Aging-associated alterations in mammary epithelia and stroma revealed by single-cell RNA sequencing. Cell Rep. 2020;33(13):108566. https://doi.org/10.1016/j.celrep.2020.108566.
Kanaya N, Chang G, Wu X, Saeki K, Bernal L, Shim HJ, et al. Single-cell RNA-sequencing analysis of estrogen- and endocrine-disrupting chemical-induced reorganization of mouse mammary gland. Commun Biol. 2019;2:406. https://doi.org/10.1038/s42003-019-0618-9.
Wu SZ, Roden DL, Wang C, Holliday H, Harvey K, Cazet AS, et al. Stromal cell diversity associated with immune evasion in human triple-negative breast cancer. Embo J. 2020;39(19):e104063. https://doi.org/10.15252/embj.2019104063.
Sun H, Miao Z, Zhang X, Chan UI, Su SM, Guo S, et al. Single-cell RNA-Seq reveals cell heterogeneity and hierarchy within mouse mammary epithelia. J Biol Chem. 2018;293(22):8315–29. https://doi.org/10.1074/jbc.RA118.002297.
McNally S, Martin F. Molecular regulators of pubertal mammary gland development. Ann Med. 2011;43(3):212–34. https://doi.org/10.3109/07853890.2011.554425.
McCave EJ, Cass CA, Burg KJ, Booth BW. The normal microenvironment directs mammary gland development. J Mammary Gland Biol Neoplasia. 2010;15(3):291–9. https://doi.org/10.1007/s10911-010-9190-0.
Kumar T, Nee K, Wei R, He S, Nguyen QH, Bai S, et al. A spatially resolved single-cell genomic atlas of the adult human breast. Nature. 2023. https://doi.org/10.1038/s41586-023-06252-9.
Schlaitz AL, Thompson J, Wong CC, Yates JR 3rd, Heald R. REEP3/4 ensure endoplasmic reticulum clearance from metaphase chromatin and proper nuclear envelope architecture. Dev Cell. 2013;26(3):315–23. https://doi.org/10.1016/j.devcel.2013.06.016.
Strzelecka PM, Ranzoni AM, Cvejic A. Dissecting human disease with single-cell omics: application in model systems and in the clinic. Dis Model Mech. 2018;11(11):dmm036525. https://doi.org/10.1242/dmm.036525.
Ren X, Zhong G, Zhang Q, Zhang L, Sun Y, Zhang Z. Reconstruction of cell spatial organization from single-cell RNA sequencing data based on ligand-receptor mediated self-assembly. Cell Res. 2020;30(9):763–78. https://doi.org/10.1038/s41422-020-0353-2.
Theocharidis G, Thomas BE, Sarkar D, Mumme HL, Pilcher WJR, Dwivedi B, et al. Single cell transcriptomic landscape of diabetic foot ulcers. Nat Commun. 2022;13:181. https://doi.org/10.1038/s41467-021-27801-8.
Wu X, Guo X, Xie C, Zhang T, Gao P, Gao T, et al. Effects of a two-meal daily feeding pattern with varied crude protein levels on growth performance and antioxidant indexes in pigs. Anim Nutr. 2016;2(4):267–70. https://doi.org/10.1016/j.aninu.2016.08.002.
Fan X, Fu Y, Zhou X, Sun L, Yang M, Wang M, et al. Single-cell transcriptome analysis reveals cell lineage specification in temporal-spatial patterns in human cortical development. Sci Adv. 2020;6(34):eaaz2978. https://doi.org/10.1126/sciadv.aaz2978.
Bai YM, Yang F, Luo P, Xie LL, Chen JH, Guan YD, et al. Single-cell transcriptomic dissection of the cellular and molecular events underlying the triclosan-induced liver fibrosis in mice. Mil Med Res. 2023;10(1):7. https://doi.org/10.1186/s40779-023-00441-3.
Wang Y, Wang R, Zhang S, Song S, Wang L. iTALK: an R package to characterize and illustrate intercellular communication. bioRxiv. 2019. https://doi.org/10.1101/507871.
Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H. SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes. Nucleic Acids Res. 2021;49(9):e50. https://doi.org/10.1093/nar/gkab043.
Mcinnes L, Healy J. UMAP: Uniform Manifold Approximation and Projection. The Journal of Open Source Software. 2018;3(29):861. https://doi.org/10.21105/joss.00861.
Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):357387.E29. https://doi.org/10.1016/j.cell.2021.04.048.
Cinti S. Pink Adipocytes. Trends Endocrinol Metab. 2018;29(9):651–66. https://doi.org/10.1016/j.tem.2018.05.007.
Fu NY, Nolan E, Lindeman GJ, Visvader JE. Stem cells and the differentiation hierarchy in mammary gland development. Physiol Rev. 2020;100(2):489–523. https://doi.org/10.1152/physrev.00040.2018.
Shackleton M, Vaillant F, Simpson KJ, Stingl J, Smyth GK, Asselin-Labat ML, et al. Generation of a functional mammary gland from a single stem cell. Nature. 2006;439(7072):84–8. https://doi.org/10.1038/nature04372.
Stingl J, Eirew P, Ricketson I, Shackleton M, Vaillant F, Choi D, et al. Purification and unique properties of mammary epithelial stem cells. Nature. 2006;439(7079):993–7. https://doi.org/10.1038/nature04496.
Witkowska-Zimny M, Kaminska-El-Hassan E. Cells of human breast milk. Cell Mol Biol Lett. 2017;22:11. https://doi.org/10.1186/s11658-017-0042-4.
Anderson BM, MacLennan MB, Hillyer LM, Ma DW. Lifelong exposure to n-3 PUFA affects pubertal mammary gland development. Appl Physiol Nutr Metab. 2014;39(6):699–706. https://doi.org/10.1139/apnm-2013-0365.
Inman JL, Robertson C, Mott JD, Bissell MJ. Mammary gland development: cell fate specification, stem cells and the microenvironment. Development. 2015;142(6):1028–42. https://doi.org/10.1242/dev.087643.
Brisken C, Ataca D. Endocrine hormones and local signals during the development of the mouse mammary gland. Wiley Interdiscip Rev Dev Biol. 2015;4(3):181–95. https://doi.org/10.1002/wdev.172.
Slepicka PF, Somasundara AVH, Dos Santos CO. The molecular basis of mammary gland development and epithelial differentiation. Semin Cell Dev Biol. 2021;114:93–112. https://doi.org/10.1016/j.semcdb.2020.09.014.
Zhao W, Shahzad K, Jiang M, Graugnard DE, Rodriguez-Zas SL, Luo J, et al. Bioinformatics and gene network analyses of the swine mammary gland transcriptome during late gestation. Bioinform Biol Insights. 2013;7:193–216. https://doi.org/10.4137/bbi.S12205.
Ji F, Hurley WL, Kim SW. Characterization of mammary gland development in pregnant gilts. J Anim Sci. 2006;84(3):579–87. https://doi.org/10.2527/2006.843579x.
Tognon CE, Sorensen PH. Targeting the insulin-like growth factor 1 receptor (IGF1R) signaling pathway for cancer therapy. Expert Opin Ther Targets. 2012;16(1):33–48. https://doi.org/10.1517/14728222.2011.638626.
Wang K, Wu P, Yang Q, Chen D, Zhou J, Jiang A, et al. Detection of selection signatures in Chinese Landrace and Yorkshire pigs based on genotyping-by-sequencing data. Front Genet. 2018;9:119. https://doi.org/10.3389/fgene.2018.00119.
Tracz AF, Szczylik C, Porta C, Czarnecka AM. Insulin-like growth factor-1 signaling in renal cell carcinoma. BMC Cancer. 2016;16:453. https://doi.org/10.1186/s12885-016-2437-4.
Muraoka-Cook RS, Feng SM, Strunk KE, Earp HS 3rd. ErbB4/HER4: role in mammary gland development, differentiation and growth inhibition. J Mammary Gland Biol Neoplasia. 2008;13(2):235–46. https://doi.org/10.1007/s10911-008-9080-x.
Vijayakumar P, Bakyaraj S, Singaravadivelan A, Vasanthakumar T, Suresh R. Meta-analysis of mammary RNA seq datasets reveals the molecular understanding of bovine lactation biology. Genome. 2019;62(7):489–501. https://doi.org/10.1139/gen-2018-0144.
Marina H, Gutiérrez-Gil B, Esteban-Blanco C, Suárez-Vega A, Pelayo R, Arranz JJ. Analysis of whole genome resequencing datasets from a worldwide sample of sheep breeds to identify potential causal mutations influencing milk composition traits. Animals (Basel). 2020;10(9):1542. https://doi.org/10.3390/ani10091542.
Lucenay KS, Doostan I, Karakas C, Bui T, Ding Z, Mills GB, et al. Cyclin E associates with the lipogenic enzyme ATP-citrate lyase to enable malignant growth of breast cancer cells. Cancer Res. 2016;76(8):2406–18. https://doi.org/10.1158/0008-5472.Can-15-1646.
Morammazi S, Masoudi AA, VaezTorshizi R, Pakdel A. Differential expression of the alpha S1 casein and beta-lactoglobulin genes in different physiological stages of the Adani goats mammary glands. Iran J Biotechnol. 2016;14(4):278–85. https://doi.org/10.15171/ijb.1171.
Sun J, Liu J, Huang B, Kan X, Chen G, Wang W, et al. Kisspeptin-10 induces β-casein synthesis via GPR54 and its downstream signaling pathways in bovine mammary epithelial cells. Int J Mol Sci. 2017;18(12):2621. https://doi.org/10.3390/ijms18122621.
Kawahara R, Niwa Y, Simizu S. Integrin β1 is an essential factor in vasculogenic mimicry of human cancer cells. Cancer Sci. 2018;109(8):2490–6. https://doi.org/10.1111/cas.13693.
Chen YY, Brown NJ, Jones R, Lewis CE, Mujamammi AH, Muthana M, et al. A peptide derived from TIMP-3 inhibits multiple angiogenic growth factor receptors and tumour growth and inflammatory arthritis in mice. Angiogenesis. 2014;17(1):207–19. https://doi.org/10.1007/s10456-013-9389-y.
Olsson AK, Dimberg A, Kreuger J, Claesson-Welsh L. VEGF receptor signalling - in control of vascular function. Nat Rev Mol Cell Biol. 2006;7(5):359–71. https://doi.org/10.1038/nrm1911.
Singh NK, Hansen DE 3rd, Kundumani-Sridharan V, Rao GN. Both Kdr and Flt1 play a vital role in hypoxia-induced Src-PLD1-PKCγ-cPLA(2) activation and retinal neovascularization. Blood. 2013;121(10):1911–23. https://doi.org/10.1182/blood-2012-03-419234.
Tammela T, Zarkada G, Wallgard E, Murtomäki A, Suchting S, Wirzenius M, et al. Blocking VEGFR-3 suppresses angiogenic sprouting and vascular network formation. Nature. 2008;454(7204):656–60. https://doi.org/10.1038/nature07083.
Dick SA, Macklin JA, Nejat S, Momen A, Clemente-Casares X, Althagafi MG, et al. Self-renewing resident cardiac macrophages limit adverse remodeling following myocardial infarction. Nat Immunol. 2019;20(1):29–39. https://doi.org/10.1038/s41590-018-0272-2.
Laakkonen P, Waltari M, Holopainen T, Takahashi T, Pytowski B, Steiner P, et al. Vascular endothelial growth factor receptor 3 is involved in tumor angiogenesis and growth. Cancer Res. 2007;67(2):593–9. https://doi.org/10.1158/0008-5472.Can-06-3567.
Anisimov A, Alitalo A, Korpisalo P, Soronen J, Kaijalainen S, Leppänen VM, et al. Activated forms of VEGF-C and VEGF-D provide improved vascular function in skeletal muscle. Circ Res. 2009;104(11):1302–12. https://doi.org/10.1161/circresaha.109.197830.
Zhang W, Zhang Y, Guo X, Zeng Z, Wu J, Liu Y, et al. Sirt1 protects endothelial cells against LPS-induced barrier dysfunction. Oxid Med Cell Longev. 2017;2017:4082102. https://doi.org/10.1155/2017/4082102.
Cheng HS, Lee JXT, Wahli W, Tan NS. Exploiting vulnerabilities of cancer by targeting nuclear receptors of stromal cells in tumor microenvironment. Mol Cancer. 2019;18:51. https://doi.org/10.1186/s12943-019-0971-9.
Nwose EU, Bwititi PT. Autophagy in diabetes pathophysiology: Oxidative damage screening as potential for therapeutic management by clinical laboratory methods. Front Cell Dev Biol. 2021;9:651776. https://doi.org/10.3389/fcell.2021.651776.
Hardbower DM, Singh K, Asim M, Verriere TG, Olivares-Villagómez D, Barry DP, et al. EGFR regulates macrophage activation and function in bacterial infection. J Clin Invest. 2016;126(9):3296–312. https://doi.org/10.1172/jci83585.
Garner JB, Chamberlain AJ, Vander Jagt C, Nguyen TTT, Mason BA, Marett LC, et al. Gene expression of the heat stress response in bovine peripheral white blood cells and milk somatic cells in vivo. Sci Rep. 2020;10:19181. https://doi.org/10.1038/s41598-020-75438-2.
Langel SN, Wark WA, Garst SN, James RE, McGilliard ML, Petersson-Wolfe CS, et al. Effect of feeding whole compared with cell-free colostrum on calf immune status: The neonatal period. J Dairy Sci. 2015;98(6):3729–40. https://doi.org/10.3168/jds.2014-8422.
Ballard O, Morrow AL. Human milk composition: nutrients and bioactive factors. Pediatr Clin North Am. 2013;60(1):49–74. https://doi.org/10.1016/j.pcl.2012.10.002.
Hassiotou F, Geddes DT. Immune cell-mediated protection of the mammary gland and the infant during breastfeeding. Adv Nutr. 2015;6(3):267–75. https://doi.org/10.3945/an.114.007377.
Zimmermann J, Macpherson AJ. Breast milk modulates transgenerational immune inheritance. Cell. 2020;181(6):1202–4. https://doi.org/10.1016/j.cell.2020.05.030.
Tseng KY, Chen YH, Lin S. Zinc finger protein ZFP36L1 promotes osteoblastic differentiation but represses adipogenic differentiation of mouse multipotent cells. Oncotarget. 2017;8(13):20588–601. https://doi.org/10.18632/oncotarget.15246.
Cen Y, Youn DY, Sauve AA. Advances in characterization of human sirtuin isoforms: chemistries, targets and therapeutic applications. Curr Med Chem. 2011;18(13):1919–35. https://doi.org/10.2174/092986711795590084.
Li D, Li F, Jiang K, Zhang M, Han R, Jiang R, et al. Integrative analysis of long noncoding RNA and mRNA reveals candidate lncRNAs responsible for meat quality at different physiological stages in Gushi chicken. PLoS One. 2019;14(4):e0215006. https://doi.org/10.1371/journal.pone.0215006.
Grimsey N, Han GS, O’Hara L, Rochford JJ, Carman GM, Siniossoglou S. Temporal and spatial regulation of the phosphatidate phosphatases lipin 1 and 2. J Biol Chem. 2008;283(43):29166–74. https://doi.org/10.1074/jbc.M804278200.
Massara M, Bonavita O, Savino B, Caronni N, MollicaPoeta V, Sironi M, et al. ACKR2 in hematopoietic precursors as a checkpoint of neutrophil release and anti-metastatic activity. Nat Commun. 2018;9:676. https://doi.org/10.1038/s41467-018-03080-8.
Ju Y, Sun C, Wang X. Loss of atypical chemokine receptor 4 facilitates C-C motif chemokine ligand 21-mediated tumor growth and invasion in nasopharyngeal carcinoma. Exp Ther Med. 2019;17(1):613–20. https://doi.org/10.3892/etm.2018.7007.
Xue D, Zheng Y, Wen J, Han J, Tuo H, Liu Y, et al. Role of chemokines in hepatocellular carcinoma (Review). Oncol Rep. 2021;45(3):809–23. https://doi.org/10.3892/or.2020.7906.
Pond AC, Bin X, Batts T, Roarty K, Hilsenbeck S, Rosen JM. Fibroblast growth factor receptor signaling is essential for normal mammary gland development and stem cell function. Stem Cells. 2013;31(1):178–89. https://doi.org/10.1002/stem.1266.
Cui E, Guo H, Shen M, Yu H, Gu D, Mao W, et al. Adiponectin inhibits migration and invasion by reversing epithelial-mesenchymal transition in non-small cell lung carcinoma. Oncol Rep. 2018;40(3):1330–8. https://doi.org/10.3892/or.2018.6523.
Lund LR, Rømer J, Thomasset N, Solberg H, Pyke C, Bissell MJ, et al. Two distinct phases of apoptosis in mammary gland involution: proteinase-independent and -dependent pathways. Development. 1996;122(1):181–93. https://doi.org/10.1242/dev.122.1.181.
Ward EP, Bartolone SN, Chancellor MB, Peters KM, Lamb LE. Proteomic analysis of bladder biopsies from interstitial cystitis/bladder pain syndrome patients with and without Hunner’s lesions reveals differences in expression of inflammatory and structural proteins. BMC Urol. 2020;20:180. https://doi.org/10.1186/s12894-020-00751-x.
Zhai S, Sun B, Zhang Y, Zhao L, Zhang L. IL-17 aggravates renal injury by promoting podocyte injury in children with primary nephrotic syndrome. Exp Ther Med. 2020;20(1):409–17. https://doi.org/10.3892/etm.2020.8698.
Hu S, Luo Q, Cun B, Hu D, Ge S, Fan X, et al. The pharmacological NF-κB inhibitor BAY11–7082 induces cell apoptosis and inhibits the migration of human uveal melanoma cells. Int J Mol Sci. 2012;13(12):15653–67. https://doi.org/10.3390/ijms131215653.
Chermuła B, Jeseta M, Sujka-Kordowska P, Konwerska A, Jankowski M, Kranc W, et al. Genes regulating hormone stimulus and response to protein signaling revealed differential expression pattern during porcine oocyte in vitro maturation, confirmed by lipid concentration. Histochem Cell Biol. 2020;154(1):77–95. https://doi.org/10.1007/s00418-020-01866-w.
Liu C, Huang W, Lei Q. Regulation and function of the TAZ transcription co-activator. Int J Biochem Mol Biol. 2011;2(3):247–56.
Perotti C, Wiedl T, Florin L, Reuter H, Moffat S, Silbermann M, et al. Characterization of mammary epithelial cell line HC11 using the NIA 15k gene array reveals potential regulators of the undifferentiated and differentiated phenotypes. Differentiation. 2009;78(5):269–82. https://doi.org/10.1016/j.diff.2009.05.003.
Lodhi IJ, Semenkovich CF. Peroxisomes: a nexus for lipid metabolism and cellular signaling. Cell Metab. 2014;19(3):380–92. https://doi.org/10.1016/j.cmet.2014.01.002.
Zwick RK, Rudolph MC, Shook BA, Holtrup B, Roth E, Lei V, et al. Adipocyte hypertrophy and lipid dynamics underlie mammary gland remodeling after lactation. Nat Commun. 2018;9:3592. https://doi.org/10.1038/s41467-018-05911-0.
de Haan W, Bhattacharjee A, Ruddle P, Kang MH, Hayden MR. ABCA1 in adipocytes regulates adipose tissue lipid content, glucose tolerance, and insulin sensitivity. J Lipid Res. 2014;55(3):516–23. https://doi.org/10.1194/jlr.M045294.
Schmitz G, Langmann T. Transcriptional regulatory networks in lipid metabolism control ABCA1 expression. Biochim Biophys Acta. 2005;1735(1):1–19. https://doi.org/10.1016/j.bbalip.2005.04.004.
Silverstein RL, Febbraio M. CD36, a scavenger receptor involved in immunity, metabolism, angiogenesis, and behavior. Sci Signal. 2009;2(72):re3. https://doi.org/10.1126/scisignal.272re3.
Wang H, Eckel RH. Lipoprotein lipase: from gene to obesity. Am J Physiol Endocrinol Metab. 2009;297(2):E271–88. https://doi.org/10.1152/ajpendo.90920.2008.
Li Y, Yu S, Chen L, Hu X, Zheng J, Deng X. Involvement of PPARγ/FSP27 in the pathogenic mechanism underlying insulin resistance: tipping the balance between lipogenesis and fat storage in adult catch-up growth rats. Nutr Metab (Lond). 2019;16:11. https://doi.org/10.1186/s12986-019-0336-9.
Arnandis T, Ferrer-Vicens I, García-Trevijano ER, Miralles VJ, García C, Torres L, et al. Calpains mediate epithelial-cell death during mammary gland involution: mitochondria and lysosomal destabilization. Cell Death Differ. 2012;19(9):1536–48. https://doi.org/10.1038/cdd.2012.46.
The authors would like to thank all members of this work for their advice and technical assistance.
This work was supported by the National Key R&D Program of China (2020YFA0509500, 2021YFD1301101 and 2021YFA0805903), the Sichuan Science and Technology Program, (2023YFN0088 and 2021YFYZ0030), the National Center of Technology Innovation for Pigs (SCCXTD-2023–08) and the National Natural Science Foundation of China (32272837 and 32225046), Tianfu Agricultural Master Project.
Ethics approval and consent to participate
All animal procedures were performed in accordance with animal management regulations and approved by the animal ethical and welfare committee (AEWC) of Sichuan Agricultural University under permit No. DKY-B20171902. All efforts were made to minimize animal suffering.
Consent for publication
The authors declare that they have no competing interests.
The phenotype data of adipocytes in the mammary gland (mean ± SD). Table S2. The information of sequencing statistics and cell statistics. Table S3. Maker genes used for cell type annotation. Table S4. Fraction of each cell type in each developmental stage. Table S5. Statistics of the number of differentially expressed genes identified in two adjacent periods in each cell type. Fig. S1. The diagrams showed the nFeatures, nCounts, mitochondrial percent and ribosomal percent of each sample before or after quality control. Fig. S2. The typical L-R interactions predicted by iTALK between any two cell types. Fig. S3. Deconvolution of ST data based on snRNA-seq data. Fig. S4. GO annotation and KEGG pathway analysis of upregulated DEGs identified in adjacent developmental stages. Fig. S5. The common significant GO terms and KEGG pathways identified in adjacent developmental stages. Fig. S6. Venn graph of the numbers of shared and unique significant GO term sets or KEGG pathway sets in different cell types (adipocytes, epithelial, fibroblasts, endothelial, myoepithelial, immune and precursor cells; indicated with different color) in the mammary gland of five devel-opmental stages (L0 vs. G90, L20 vs. L0, PI2 vs. L20, and PI7 vs. PI2 groups).
About this article
Cite this article
Fan, Y., Jin, L., He, Z. et al. A cell transcriptomic profile provides insights into adipocytes of porcine mammary gland across development. J Animal Sci Biotechnol 14, 126 (2023). https://doi.org/10.1186/s40104-023-00926-0