Advertisement

In vitro transdifferentiated signatures of goat preadipocytes into mammary epithelial cells revealed by DNA methylation and transcriptome profiling

Open AccessPublished:October 15, 2022DOI:https://doi.org/10.1016/j.jbc.2022.102604
      During mammary development, the transdifferentiation of mammary preadipocytes is one of the important sources for lactating mammary epithelial cells (MECs). However, there is limited knowledge about the mechanisms of dynamic regulation of transcriptome and genome-wide DNA methylation in the preadipocyte transdifferentiation process. Here, to gain more insight into these mechanisms, preadipocytes were isolated from adipose tissues from around the goat mammary gland (GM-preadipocytes). The GM-preadipocytes were cultured on Matrigel in conditioned media made from goat MECs to induce GM-preadipocyte-to-MEC transdifferentiation. The transdifferentiated GM-preadipocytes showed high abundance of keratin 18, which is a marker protein of MECs, and formed mammary acinar-like structures after 8 days of induction. Then, we performed transcriptome and DNA methylome profiling of the GM-preadipocytes and transdifferentiated GM-preadipocytes, respectively, and the differentially expressed genes and differentially methylated genes that play underlying roles in the process of transdifferentiation were obtained. Subsequently, we identified the candidate transcription factors in regulating the GM-preadipocyte-to-MEC transdifferentiation by transcription factor–binding motif enrichment analysis of differentially expressed genes and differentially methylated genes. Meanwhile, the secretory proteome of GM-preadipocytes cultured in conditioned media was also detected. By integrating the transcriptome, DNA methylome, and proteome, three candidate genes, four proteins, and several epigenetic regulatory axes were further identified, which are involved in regulation of the cell cycle, cell polarity establishment, cell adhesion, cell reprogramming, and adipocyte plasticity. These findings provide novel insights into the molecular mechanism of preadipocyte transdifferentiation and mammary development.

      Keywords

      Abbreviations:

      ACTG1 (actin gamma 1), ALB (albumin), ASC (adipose-derived stem cell), BACH2 (BTB domain and CNC homolog 2), cDNA (complementary DNA), CM (conditioned media), CTCF (CCCTC-binding factor), DEG (differentially expressed gene), DLK1 (delta-like noncanonical Notch ligand 1), DMEM (Dulbecco's modified Eagle's medium), DMG (differentially methylated gene), DML (differentially methylated location), DMR (differentially methylated region), ECM (extracellular matrix), FBS (fetal bovine serum), GMEC (goat MEC), GO (Gene Ontology), GM-preadipocyte (goat mammary preadipocyte), HSP90AA1 (heat shock protein 90 alpha family class A member 1), KEGG (Kyoto Encyclopedia of Genes and Genomes), KRT18 (keratin 18), MEC (mammary epithelial cell), NG (gene that has significantly negative correlation), NID1 (nidogen 1), PDGFRA+ (platelet-derived growth factor receptor alpha (+)), PLCL1 (phospholipase C like 1), PPAR (peroxisome proliferator–activated receptor), RND1 (Rho family GTPase 1), RT–qPCR (real-time quantitative PCR), SPPL3 (signal peptide peptidase like 3), STAT (signal transducer and activator of transcription), TEM (transmission electron microscope), TF (transcription factor), WGBS (whole genome bisulfite sequencing)
      The mammary gland is a tissue that develops throughout the sexual maturity and performs physiological function of milk synthesis and secretion. It is mostly made up of epithelial tissue and stroma, which contains adipose tissue known as the “fat pad.” The fat pad in mammary gland is required for mammary development via signals that trigger the ductal morphogenesis (
      • Neville M.C.
      • Medina D.
      • Monks J.
      • Hovey R.C.
      The mammary fat pad.
      ). There are three kinds of adipocytes in fat pad, white adipocytes that store and breakdown lipids, brown adipocytes that are thermogenic, and pink adipocytes with lactation function (
      • Giordano A.
      • Perugini J.
      • Kristensen D.M.
      • Sartini L.
      • Frontini A.
      • Kajimura S.
      • et al.
      Mammary alveolar epithelial cells convert to brown adipocytes in post-lactating mice.
      ). Under certain circumstances, all three types of adipocytes have been shown to coexpress leptin and perilipin A and transform accordingly (
      • Giordano A.
      • Smorlesi A.
      • Frontini A.
      • Barbatelli G.
      • Cinti S.
      White, brown and pink adipocytes: the extraordinary plasticity of the adipose organ.
      ). When hormone levels fluctuate during different stages of development, they may also go through a series of modifications in the mammary gland (
      • Smith G.H.
      • Medina D.
      Re-evaluation of mammary stem cell biology based on in vivo transplantation.
      ). Furthermore, adipocytes remain active throughout mammary development, demonstrating a significant amount of flexibility while transitioning from adipocytes to fibroblast-like cells and preadipocytes. The mRNA and protein levels of delta-like noncanonical Notch ligand 1 (DLK1) are highly expressed in preadipocytes and absent in mature adipocytes (
      • Joshi P.A.
      • Waterhouse P.D.
      • Kasaian K.
      • Fang H.
      • Gulyaeva O.
      • Sul H.S.
      • et al.
      PDGFRalpha(+) stromal adipocyte progenitors transition into epithelial cells during lobulo-alveologenesis in the murine mammary gland.
      ,
      • Colleluori G.
      • Perugini J.
      • Barbatelli G.
      • Cinti S.
      Mammary gland adipocytes in lactation cycle, obesity and breast cancer.
      ,
      • Wang Q.A.
      • Song A.
      • Chen W.
      • Schwalie P.C.
      • Zhang F.
      • Vishvanath L.
      • et al.
      Reversible de-differentiation of mature white adipocytes into preadipocyte-like precursors during lactation.
      ). The mRNA and protein levels of CD34 molecule are highly expressed in preadipocytes and lowly expressed in mature adipocytes (
      • Lee Y.H.
      • Petkova A.P.
      • Mottillo E.P.
      • Granneman J.G.
      In vivo identification of bipotential adipocyte progenitors recruited by beta3-adrenoceptor activation and high-fat feeding.
      ,
      • Carriere A.
      • Jeanson Y.
      • Cote J.A.
      • Dromard C.
      • Galinier A.
      • Menzel S.
      • et al.
      Identification of the ectoenzyme CD38 as a marker of committed preadipocytes.
      ). Hence, DLK1 and CD34 have been utilized as marker proteins of preadipocytes. The human preadipocytes marked with DLK1 cultured in the medium made from human mammary epithelial cells (MECs) in 3D culture condition transdifferentiated into adult MECs (
      • Joshi P.A.
      • Waterhouse P.D.
      • Kasaian K.
      • Fang H.
      • Gulyaeva O.
      • Sul H.S.
      • et al.
      PDGFRalpha(+) stromal adipocyte progenitors transition into epithelial cells during lobulo-alveologenesis in the murine mammary gland.
      ). Preadipocytes tagged with platelet-derived growth factor receptor alpha (+) (PDGFRA+), as a source of epithelial descendents, depleted during mammary acini production after transdifferentiating into MECs (also known as “pink adipocytes”) (
      • Joshi P.A.
      • Waterhouse P.D.
      • Kasaian K.
      • Fang H.
      • Gulyaeva O.
      • Sul H.S.
      • et al.
      PDGFRalpha(+) stromal adipocyte progenitors transition into epithelial cells during lobulo-alveologenesis in the murine mammary gland.
      ,
      • Colleluori G.
      • Perugini J.
      • Barbatelli G.
      • Cinti S.
      Mammary gland adipocytes in lactation cycle, obesity and breast cancer.
      ,
      • Wang Q.A.
      • Song A.
      • Chen W.
      • Schwalie P.C.
      • Zhang F.
      • Vishvanath L.
      • et al.
      Reversible de-differentiation of mature white adipocytes into preadipocyte-like precursors during lactation.
      ,
      • Berry R.
      • Rodeheffer M.S.
      Characterization of the adipocyte cellular lineage in vivo.
      ). Keratin 18 (KRT18) is highly expressed in monolayer or simple epithelial tissue and was used as a marker protein of MECs (
      • Faridi N.
      • Bathaie S.Z.
      • Abroun S.
      • Farzaneh P.
      • Karbasian H.
      • Tamanoi F.
      • et al.
      Isolation and characterization of the primary epithelial breast cancer cells and the adjacent normal epithelial cells from Iranian women's breast cancer tumors.
      ,
      • Weng Y.R.
      • Cui Y.
      • Fang J.Y.
      Biological functions of cytokeratin 18 in cancer.
      ).
      Acinus is a spherical body of monolayer MECs that surrounds center lumen, and these acini in the mammary gland are involved in lactation. The production of acinus occurs in two ways: (i) the proliferation and differentiation of mammary epithelial stem cells, which are associated with expansion and maturation of the epithelial portion and (ii) pregnancy-induced transdifferentiation of mammary adipocytes to MECs (
      • Hovey R.C.
      • Trott J.F.
      Morphogenesis of mammary gland development.
      ). MECs are derived from acini progenitor cells in the first stage of pregnancy, and lipid-laden MECs develop in the second stage when subcutaneous adipocytes shrink and exhibit epithelial-like characteristics to transdifferentiate into MECs. However, these cells produce the secretory acini with myoepithelial cells (
      • Giordano A.
      • Smorlesi A.
      • Frontini A.
      • Barbatelli G.
      • Cinti S.
      White, brown and pink adipocytes: the extraordinary plasticity of the adipose organ.
      ,
      • Yin Y.
      • Yuan H.
      • Wang C.
      • Pattabiraman N.
      • Rao M.
      • Pestell R.G.
      • et al.
      3-phosphoinositide-dependent protein kinase-1 activates the peroxisome proliferator-activated receptor-gamma and promotes adipocyte differentiation.
      ) that mature in three consecutive steps. The first is proliferation, in which single cells grow into multicellular clusters; the second is secession, in which outer layer of cells establishes a polarity axis and undergoes proliferative suppression; and the third is apoptosis, in which inner cells die, resulting in the formation of hollow lumen (
      • Whyte J.
      • Thornton L.
      • McNally S.
      • McCarthy S.
      • Lanigan F.
      • Gallagher W.M.
      • et al.
      PKCzeta regulates cell polarisation and proliferation restriction during mammary acinus formation.
      • Muthuswamy S.K.
      • Li D.
      • Lelievre S.
      • Bissell M.J.
      • Brugge J.S.
      ErbB2, but not ErbB1, reinitiates proliferation and induces luminal repopulation in epithelial acini.
      ). The transdifferentiated preadipocytes also undergo the process of acinar formation. We have known that preadipocyte-to-MEC transdifferentiation is crucial for mammary development because it also involves in the production of acinar structure. However, earlier work has revealed little about the processes of dynamic control of transcriptome and genome-wide DNA methylation in preadipocyte transdifferentiation, and there has been little in-depth investigation into the molecular mechanisms involved in preadipocyte transdifferentiation.
      Although mouse models of mammary development are very valuable and easy to handle, they have inherent limitations, and caution is recommended while applying the results of development studies directly to humans (
      • Gusterson B.A.
      • Stein T.
      Human breast development.
      ). Previous research has indicated that ruminant mammary could be utilized as a good supplementary model for the study of mammary development because of in-depth understanding development process and microanatomical of ruminant mammary (
      • Nagy D.
      • Gillis C.M.C.
      • Davies K.
      • Fowden A.L.
      • Rees P.
      • Wills J.W.
      • et al.
      Developing ovine mammary terminal duct lobular units have a dynamic mucosal and stromal immune microenvironment.
      ). Therefore, in this study, the preadipocytes were isolated from adipose tissues from around the goat mammary gland (GM-preadipocytes), and the cultured medium of goat MECs (GMECs) grown on Matrigel was collected as the conditioned media (CM) to induce GM-preadipocyte-to-MEC transdifferentiation. Furthermore, GM-preadipocytes and transdifferentiated GM-preadipocytes performed transcriptome and DNA methylome sequencing to screen differentially expressed genes (DEGs) and differentially methylated genes (DMGs). Moreover, the proteome of CM was also detected. The candidate genes, proteins, and epigenetic regulatory networks involved in GM-preadipocyte-to-MEC transdifferentiation were further identified through integration analysis of the transcriptome, DNA methylome, and the proteome. This work will provide novel insights into molecular mechanism of adipocyte plasticity and supply a theoretical underpinning for the mammary development.

      Results

      The characterization of GMECs and GM-preadipocytes

      The GM-preadipocytes were distinguished from GMECs by Oil Red O staining and immunofluorescence of marker proteins. Preadipocyte marker proteins included DLK1 and CD34, whereas epithelial cell marker protein included KRT18. During adipogenic induction, the growth medium of GM-preadipocytes was changed to adipogenic induction medium at day 0. A small number of lipid droplets in cells were observed at day 2. The GM-preadipocytes have differentiation potential, as indicated by the increase in lipid droplet formation at day 5 and the remarkable increase in droplet size that continued until day 8 (Fig. S1A). Moreover, CD34 and DLK1 were highly expressed in GM-preadipocytes, and KRT18 was found to be expressed strongly in GMECs (Fig. S1B). These findings showed that GM-preadipocytes and GMECs may be utilized for future researches.

      Morphological changes of the GMECs cultured in Matrigel

      GMECs were cultured with the growth medium in Matrigel to simulate the conditions seen in vivo to determine whether they might spontaneously induce to the formation of acinar structures. During spontaneous induction of GMECs, the cell suspension was inoculated on Matrigel at day 0. Then we observed cells regularly and stably embedded in Matrigel at day 2; and cell clusters gradually formed at day 4, suggesting the establishment of polarity of outer cells and reduction of proliferative activity; at day 6, cell clusters got bigger, demonstrating the presence of cell adhesion; and at day 8, acinar structures and lumen appeared, indicating that the inner cells suffered apoptosis (Fig. 1). These findings showed that the GMECs cultured on Matrigel-coated dishes produced acinar structures and lumens.
      Figure thumbnail gr1
      Figure 1The morphological characteristics of the epithelial structures formed by GMECs cultured in Matrigel. The morphological changes of GMECs at day 2, 4, 6, and 8 of spontaneous induction, and the cell suspension was inoculated on Matrigel at day 0. Cells were embedded in Matrigel at day 2, and at day 4, cell clusters were gradually formed. At day 6, cell clusters got bigger, and at day 8, acinar structures and lumen appeared. The red boxes represent three acinar structures, respectively (magnification: 4×), and the red arrows show magnifications of the separate acinar structure (magnification: 10×). The orange and purple boxes represent different parts of an acinar structure, respectively (magnification: 10×), and the arrows with corresponding colors show magnifications of the upper and lower parts of the acinar structure (magnification: 20×), magnification in brackets. GMEC, goat mammary epithelial cell.

      Different behaviors of the GM-preadipocytes induced by the CM

      The CM were utilized to induce GM-preadipocyte-to-MEC transdifferentiation in Matrigel to determine the effect of secretory functional factors from GMEC culture medium on GM-preadipocytes. The cell morphology was not significantly changed in Preadipocyte group (Fig. 2A). However, in the Gpreadipocyte group, when cells were induced with CM, cells were reticular at day 0; the antennae of cells shrank at day 4, indicating epithelialization of GM-preadipocytes. At day 8, cells formed colonies, acinar structures, and lumens in a swirling pattern, indicating that the GM-preadipocytes experienced polarity establishment, proliferative suppression, cell adhesion, and apoptosis of inner cells (Fig. 2B). Interestingly, the acinar structures generated from GM-preadipocytes in Gpreadipocyte group were comparable to those of GMECs.
      Figure thumbnail gr2
      Figure 2Different behaviors of GM-preadipocyte induced by the CM. A, morphological changes of GM-preadipocyte in the Preadipocyte group (magnification: 10×). B, morphological changes of GM-preadipocyte in the Gpreadipocyte group. The cells were reticular at day 0; the antennae of cells shrank at day 4; the cells were developed in a swirling pattern to form colonies, acinar structure, and lumen at day 8 (magnification: 10×). The red and green boxes represent the cell clusters in the process of GM-preadipocyte transdifferentiation, respectively (magnification: 10×), and the arrows with corresponding colors show magnifications of the indicated red and green rectangular region (magnification: 20×), magnification in brackets. Day 0 refers to the day when GM-preadipocytes undergone serum starvation and were treated with CM. C, immunofluorescence of Gpreadipocyte group using KRT18 monoclonal antibody (green) under confocal laser scanning microscope. D, the subcellular structures of the GMEC, Preadipocyte, and Gpreadipocyte groups. The yellow and red arrows represent desmosome joining adjacent cells and dense apical granules, respectively. CM, conditioned media; er, endoplasmic reticulum; GMEC, goat mammary epithelial cell; GM-preadipocyte, goat mammary preadipocyte; KRT18, keratin 18; m, mitochondrion; mv, microvilli; n, nucleus.
      Simultaneously, immunofluorescence results revealed that GM-preadipocytes transdifferentiated into MECs with high KRT18 expression (Fig. 2C). Transmission electron microscope (TEM) was used to observe if transdifferentiated GM-preadipocytes established a polarity axis and formed acinar structures. In the Gpreadipocyte group, there were many microvilli, acrosome particles, and desmosomes that connected adjacent cells in the cytoplasm, and a significant number of mitochondria, which were identical to GMEC organelles (Fig. 2D). These findings showed that CM enabled GM-preadipocytes to transdifferentiate into MECs, establish polarity axis, and form acinar structures.

      Transcriptome profiling of the GM-preadipocytes pre- and post-transdifferentiation

      The transcriptome sequencing of GM-preadipocytes and transdifferentiated GM-preadipocytes was used to screen candidate genes that involved in the transdifferentiation. The sequencing results showed an average of 44,475,402 and 42,430,324 raw reads generated, respectively, whereas 42,691,545 and 41,727,560 clean reads were obtained accordingly after removing low-quality reads, with about 90.85% and 90.63% of clean reads independently mapped to goat reference genome ARS1 in the Preadipocyte and Gpreadipocyte groups, respectively (Table 1).
      Table 1The summary of data generated by transcriptome sequencing
      Sample IDRaw data (G)Raw readsClean readsQ30 (%)Mapping rate (%)
      Control 16.6944,622,18642,663,98092.8189.76
      Control 26.2141,431,86839,394,43293.2291.06
      Control 36.9846,511,60044,716,41293.1291.14
      Control 46.6644,409,91843,179,25293.3191.23
      Control 56.8145,401,43643,503,64893.3791.04
      Treatment 16.6144,034,26643,399,38693.2689.94
      Treatment 26.2341,529,49240,715,22693.4092.05
      Treatment 35.9839,863,25639,005,84093.0489.50
      Treatment 46.2241,488,54440,888,23492.8591.50
      Treatment 56.7945,236,06444,629,11692.4290.14
      Control: the Preadipocyte group; treatment: the Gpreadipocyte group.
      There were no outlier samples in either of two groups according to the principal component analysis (Fig. 3A). The PDGFRA expression found in the transcriptome of Preadipocyte group indicated that the GM-preadipocytes were PDGFRA (+) preadipocytes. The mRNA profiles of Preadipocyte and Gpreadipocyte groups indicated 823 and 887 stage-specific genes, respectively (Table S1). In comparison to the Preadipocyte group, the Gpreadipocyte group had 4337 DEGs, of which 2238 were upregulated, including several epithelial marker genes, such as factor V/VIII domain containing (MFGE8), perilipin 2 (Plin2), and E74 like ETS transcription factor (TF) 5 (Elf5); and 2099 were downregulated, including adipocyte marker genes like peroxisome proliferator–activated receptor gamma (PPARγ) and FABP4 (Fig. 3B and Table S2). Subsequently, the 10 genes selected from all DEGs were verified via real-time quantitative PCR (RT–qPCR). The results were in accordance with the transcriptome data, indicating that transcriptome sequencing generated reliable data (Fig. S2).
      Figure thumbnail gr3
      Figure 3Transcriptomic profiling of the GM-preadipocytes pre- and post-transdifferentiation. A, PCA of transcriptome. B, volcano plots of mRNA. C, KEGG pathway enrichment analysis of stage-specific genes. Venn diagram represents the intersection of pathways enriched by stage-specific genes in control and treatment groups. Control and treatment groups represent the Preadipocyte and Gpreadipocyte groups, respectively. GM-preadipocyte, goat mammary preadipocyte; KEGG, Kyoto Encyclopedia of Genes and Genomes; PCA, principal component analysis.
      Furthermore, the functions of stage-specific genes were investigated by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, and the stage-specific genes of the Preadipocyte group were enriched in PPAR, Jak- signal transducer and activator of transcription (STAT), Wnt signaling pathways, and the regulation of pluripotency of stem cell pathway. Whereas, the stage-specific genes of the Gpreadipocyte group were enriched in tight junction, gap junction, Rap1, MAPK, PI3K-Akt, and Hippo signaling pathways, which highlighted the central roles of the cell adhesion, morphogenesis, and organogenesis during the process of transdifferentiation (Fig. 3C and Table S3). Among them, both PPAR and Hippo signaling pathways have been proved to associate with the adipocyte transdifferentiation and acinar formation.
      Simultaneously, gene set enrichment analysis was used to evaluate the potential function of DEGs, and results showed that downregulated and upregulated DEGs were enriched in different Gene Ontology (GO) terms. The downregulated DEGs were significantly enriched in GO terms including focal adhesion, actin cytoskeleton organization, DNA-binding TF binding, and transcription coactivator activity. Whereas upregulated DEGs were evidently enriched in GO terms including catalytic activity, endosome, and Golgi apparatus (Fig. 4A and Table S4). Moreover, results from KEGG pathway enrichment analysis showed that upregulated DEGs were significantly involved in cell cycle and cell adhesion, whereas downregulated DEGs were strongly associated with the immune response, cell cycle, cell adhesion, morphogenesis, and organogenesis (Fig. 4B and Table S5). The overall GO and KEGG enrichment analysis demonstrated that identified DEGs were closely related to regulation of cell cycle, cell adhesion, and morphogenesis.
      Figure thumbnail gr4
      Figure 4Enrichment analysis of GSEA, KEGG, and TF motifs for DEGs. A, GSEA of DEGs. B, KEGG pathway enrichment analysis of DEGs. C and D, TF-binding motif enrichment analysis around the promoter of downregulated and upregulated DEGs, respectively. DEG, differentially expressed gene; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; TF, transcription factor.
      The promoter motif enrichment analysis was used to determine if DEGs were collectively controlled by certain TFs. We found 11 and 13 binding motifs for downregulated and upregulated DEGs, respectively. These discovered TFs might play critical roles in the various processes of GM-preadipocyte-to-MEC transdifferentiation (Fig. 4, C, D and Table S6).

      DNA methylation patterns of the GM-preadipocytes pre- and post-transdifferentiation

      The DNA methylation patterns in the GM-preadipocyte-to-MEC transdifferentiation were studied using whole genome bisulfite sequencing (WGBS). In the Preadipocyte and Gpreadipocyte groups, an average of 328,480,510 and 349,232,847 raw reads were generated, 227,236,312 and 225,855,304 clean reads were obtained, about 93.38% and 93.21% of clean reads were uniquely mapped to goat reference genome ARS1, respectively (Table 2).
      Table 2The summary of data generated by WGBS
      Sample IDRaw data (G)Raw readsClean readsMapping rate (%)
      Control 167.36336,798,123244,789,29493.22
      Control 261.26306,348,087184,612,12193.42
      Control 370.12350,612,933252,307,52093.50
      Treatment 162.42312,090,750198,212,07293.22
      Treatment 285.96429,888,384301,301,90693.21
      Treatment 361.14305,719,408178,051,93493.21
      Control: the Preadipocyte group; treatment: the Gpreadipocyte group.
      Besides, principal component analysis results revealed no outliers in all samples (Fig. S3). Moreover, the methylation degrees of different sites on each chromosome were obtained (Table S7). Average CpG counts were 7,639,657 and 5,039,938, and methylated CpG represented 63.97% and 62.83% in the Preadipocyte and Gpreadipocyte groups, respectively (Table 3). A total of 577 differentially methylated locations (DMLs) and 17 differentially methylated regions (DMRs) were identified and annotated, yielding a total of 152 DMGs (Fig. 5A and Table S8).
      Table 3Information about methylation
      Sample IDMethylated CpGTotal CpGMethylated/total CpG (%)Methylated CHTotal CHMethylated/total CH (%)
      Control 16,492,6219,731,99466.712,135,810259,575,2460.82
      Control 23,546,8625,684,47662.401,491,918162,976,2950.92
      Control 34,712,4837,502,50062.811,801,088224,005,1060.80
      Treatment 13,066,1814,861,05563.081,214,285146,354,7230.83
      Treatment 23,160,5075,087,14962.131,394,183154,033,4000.91
      Treatment 33,272,4675,171,61163.281,096,083138,180,1710.79
      Control: the Preadipocyte group; treatment: the Gpreadipocyte group.
      Figure thumbnail gr5
      Figure 5DNA methylation profiling of the GM-preadipocytes pre- and post-transdifferentiation. A, Manhattan map of DMGs, including the methylation degree of DMLs and the number of methylation sites in DMRs on every chromosome. B, GO and KEGG enrichment analysis of DMGs. C, TF-binding motif enrichment analysis around the promoter of DMGs. D, the DNA methylation–meditated TF-gene regulatory network of GM-preadipocyte-to-MEC transdifferentiation. DMG, differentially methylated gene; DML, differentially methylated location; DMR, differentially methylated region; GM-preadipocyte, goat mammary preadipocyte; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MEC, mammary epithelial cell; TF, transcription factor.
      Both GO and KEGG enrichment analyses were performed to investigate the potential functions of DMGs. The DMGs were mostly enriched in the GO terms of actin cytoskeleton, adherens junction, protein homodimerization activity, and signaling receptor binding, and KEGG pathways including phospholipase D signaling pathway, lysine degradation pathway, cell adhesion molecules, PI3K–Akt, and Jak–STAT signaling pathway (Fig. 5B and Table S9). These results confirmed that there are significant differences between the GM-preadipocytes and transdifferentiated GM-preadipocytes in both the transcriptome and DNA methylome levels. The function of DMGs was similar to that of DEGs, mainly focusing on cell cycle, cell adhesion, maintenance of cell architecture, and stability.
      The gene expression was affected with DNA methylation by modulating bindings of upstream TFs. In this study, a total of 25 DMGs that were methylated in promoter region were screened, and these promoters were highly enriched in binding motifs of different TFs, including BTB domain and CNC homolog 2 (BACH2), CCCTC-binding factor (CTCF), and Spi-1 proto-oncogene (Fig. 5C and Table S10). The interactions of TF with their target genes were then used to build the TF-gene regulatory networks. We identified links between 14 TFs and 15 DMGs and constructed a putative regulatory network associated with GM-preadipocyte transdifferentiation (Fig. 5D).

      DNA methylation regulates GM-preadipocyte transdifferentiation by affecting gene expression

      Pearson correlation coefficients were calculated to assess the relationship between gene expression and DNA methylation level of corresponding genes during the GM-preadipocyte-to-MEC transdifferentiation (Table S11). A total of 37 genes had significantly positive correlation (PGs), of which were related to cell adhesion, cell–cell junction, cell migration, and regulation of cell shape. We further found that 32 genes showed a significantly negative correlation (NGs) between DNA methylation level and gene expression, which were associated with planar polarity, regulation of cell adhesion, mammary gland epithelial cell differentiation, mammary gland epithelium development, and mammary gland alveolus development (Fig. 6A and Table S12). It seemed that function of NGs was closely related to the GM-preadipocyte-to-MEC transdifferentiation. Following that, through integration of the DEGs and DMGs, we found that 29 genes were methylated resulting in differential expression, and these genes were annotated into lipid metabolism, cell adhesion, and morphogenesis (Fig. 6B and Table S13).
      Figure thumbnail gr6
      Figure 6DNA methylation regulates GM-preadipocyte transdifferentiation by affecting gene expression. A, the distribution of correlation coefficients between DNA methylation and the expression level of corresponding genes, and GO annotation of PGs or NGs. B, intersection of DEGs and DMGs, and their GO annotation. The outer to inner circles represent chromosome names, DEGs, and DMGs, respectively. C, the DNA methylation site distribution of three candidate DMGs. The differentially DNA methylated sites are indicated by red number. The bar chart represents the comparison of DNA methylation levels between two groups for each site. Control and treatment groups represent the Preadipocyte and Gpreadipocyte group, respectively. DEG, differentially expressed gene; DMG, differentially methylated gene; GM-preadipocyte, goat mammary preadipocyte; GO, Gene Ontology; NG, gene that has significantly negative correlation; PG, genes that have significantly positive correlation.
      Based on correlation analysis and integration of DEGs with DMGs, three DMGs were identified as being involved in methylation-mediated transdifferentiation: signal peptide peptidase like 3 (SPPL3), Rho family GTPase 1 (RND1), and phospholipase C like 1 (PLCL1). The DNA methylation site distribution of these three genes showed that differential methylation occurred at eight sites in the promoter of SPPL3, and the site of 8,004,918 location was most significant, whereas the site of 30,564,254 location had significant differential methylation in the promoter of RND1, and for promoter of PLCL1, the site of 49,627,299 location was also significantly differentially methylated (Fig. 6C). These differentially methylated sites may be important for epigenetic changes of GM-preadipocyte-to-MEC transdifferentiation.

      Construction of epigenetic regulatory networks

      The GM-preadipocytes proceed through epithelization, polarity establishment, and acinar development in the presence of CM. Protein mass spectrometry was used to detect proteins in CM to explore functional factors, and 32 proteins were identified (Fig. 7A and Table S14). Functions of these proteins were annotated into different GO terms (Fig. 7B and Table S15). Among the terms, the establishing cell polarity and morphogenesis of a polarized epithelium showed that secreted proteins could regulate the polarity of GM-preadipocytes. The cell adhesion, cell–cell junctions, and positively regulated cell substrate adhesion were associated with cell cycle and acinar formation. For identified proteins, albumin (ALB), nidogen 1 (NID1), actin gamma 1 (ACTG1), and heat shock protein 90 alpha family class A member 1 (HSP90AA1) were annotated into regulation of apoptosis process, cellular process, morphogenesis of a polarized epithelium, and establishment of cell polarity, respectively. Therefore, the results suggested that ALB, NID1, ACTG1, and HSP90AA1 are essential for regulation of different transdifferentiation processes.
      Figure thumbnail gr7
      Figure 7The epigenetic regulation network for inducing GM-preadipocyte-to-MEC transdifferentiation. A, flowchart of secretory proteome detection. B, functional annotation of secretory protein in the CM. C, the construction of epigenetic regulation networks. The green, blue, and pink circles represent the DNA methylase, secretory protein of the CM, and DMGs, respectively. D, regulatory axis related to the GM-preadipocyte-to-MEC transdifferentiation. The connection in the middle is the corresponding relationship. CM, conditioned media; DMG, differentially methylated gene; GM, goat mammary; MEC, mammary epithelial cell.
      The spectrum of transdifferentiation factors was then limited to paracrine mediators, which transfer signals from GMECs to GM-preadipocytes to induce transdifferentiation. An epigenetic regulatory network was constructed by integrating the DNA methylase, DMGs, and secretory proteins (Fig. 7C and Table S16). In the network, ALB was secreted from GMECs and impacted genomic methylation by interacting directly with DNA methylase and DMGs. In addition, ALB regulated epithelial cell proliferation, apoptosis, and autophagy by interacting with cyclin-dependent kinase 6 (CDK6) and activating TF 6 (ATF6), whereas ACTG1 regulated cell junction, cell shape, cell process, and cell activity by interacting with WD repeat domain 1 (WDR1), abl interactor 1 (ABI1), and myosin VA (MYO5A). HSP90AA1 regulated cell polarity establishment and cell junction by interacting with ABI1, and NID1 regulated cell adhesion by interacting with RND1 (Fig. 7D).
      Based on four kinds of analysis, including the correlation analysis, the overlap of DEGs and DMGs, TF regulatory network of DMGs, and epigenetic regulatory network, we found that TFs CTCF and BACH2 were the central regulators and involved in the transdifferentiation of GM-preadipocytes. For example, NID1 affected binding of CTCF to RND1 by regulating the DNA methylation of RND1 and ACTG1 affects the binding of BACH2 and MYO5A by mediating the DNA methylation of MYO5A.

      Discussion

      The lack of polarized shape in secreted MECs may result in nonfunctional differentiation after disconnecting from the natural milieu. Tissue-specific genes can only be seen in vitro when extracellular matrix (ECM) signals can be received. In 3D culture, cell morphology, polarity, signal transduction, gene expression, and metabolism are similar to those of cells in vivo. The 3D cultured system provides an appropriate physiological condition for the study of complex cell–cell and cell–ECM interactions, which has been widely applied in various researches regarding mammary biology. In 3D culture condition, mouse MECs were used to study catheter invasion and elongation, morphogenetic procedures of alveolar genesis, and functional differentiation (
      • Lee E.Y.
      • Lee W.H.
      • Kaetzel C.S.
      • Parry G.
      • Bissell M.J.
      Interaction of mouse mammary epithelial cells with collagen substrata: regulation of casein gene expression and secretion.
      ,
      • Sumbal J.
      • Chiche A.
      • Charifou E.
      • Koledova Z.
      • Li H.
      Primary mammary organoid model of lactation and involution.
      ), and supplement of key components induced milk synthesis of bovine MECs (
      • Mackenzie D.D.
      • Forsyth I.A.
      • Brooker B.E.
      • Turvey A.
      Culture of bovine mammary epithelial cells on collagen gels.
      ,
      • Riley L.G.
      • Gardiner-Garden M.
      • Thomson P.C.
      • Wynn P.C.
      • Williamson P.
      • Raadsma H.W.
      • et al.
      The influence of extracellular matrix and prolactin on global gene expression profiles of primary bovine mammary epithelial cells in vitro.
      ). BME-UV1 (bovine MECs) formed acinar structures composed of polarized MECs when cells grew in Matrigel for 16 days (
      • Kozlowski M.
      • Gajewska M.
      • Majewska A.
      • Jank M.
      • Motyl T.
      Differences in growth and transcriptomic profile of bovine mammary epithelial monolayer and three-dimensional cell cultures.
      ). Many previous studies have indicated that Matrigel is an efficient 3D culture system for MECs of mouse, bovine, and human. In this study, the GMEC suspension was inoculated on Matrigel at day 0. During spontaneous induction of GMECs, we observed that cells were regularly and stably embedded in Matrigel at day 2, and the observation of cell morphology was continued every day after that in live cell imaging system. At day 4, cells rotated to form clusters and kept rotating and formed larger cell clusters at day 6. Up to day 8, many acinar structures could be observed, whereas we continuously cultured the cells for 14 days, morphology remained unchanged from day 8 to day 14. The GMEC morphology at these different time points in our results is similar to those of previous studies (
      • Wang H.
      • Lacoche S.
      • Huang L.
      • Xue B.
      • Muthuswamy S.K.
      Rotational motion during three-dimensional morphogenesis of mammary epithelial acini relates to laminin matrix assembly.
      ). Acinar structure efficiently formed from GMECs cultured in Matrigel in vitro was the experimental basis for the GM-preadipocyte-to-MEC transdifferentiation.
      This study further provided a new understanding of adipocyte plasticity, revealing tissue's complexity. Browning is an example of adipose tissue adaptation in which brown adipocytes formed in a white adipose tissue (
      • Cinti S.
      Transdifferentiation properties of adipocytes in the adipose organ.
      ,
      • Seale P.
      • Conroe H.M.
      • Estall J.
      • Kajimura S.
      • Frontini A.
      • Ishibashi J.
      • et al.
      Prdm16 determines the thermogenic program of subcutaneous white adipose tissue in mice.
      ). The formation of fat pad is required for development of mammary epithelium, which presents another unique option of adipocyte plasticity. During the pregnancy and lactation periods, a substantial number of MECs with lactation function were grown in the mammary gland, out of which 30% were derived from ductal stem cells and 70% were from transdifferentiation of adipocytes, and indicate that mammary adipocytes are able to transdifferentiate into MECs with high expression of epithelial-specific protein (
      • Morroni M.
      • Giordano A.
      • Zingaretti M.C.
      • Boiani R.
      • De Matteis R.
      • Kahn B.B.
      • et al.
      Reversible transdifferentiation of secretory epithelial cells into adipocytes in the mammary gland.
      ,
      • De Matteis R.
      • Zingaretti M.C.
      • Murano I.
      • Vitali A.
      • Frontini A.
      • Giannulis I.
      • et al.
      In vivo physiological transdifferentiation of adult adipose cells.
      ,
      • Prokesch A.
      • Smorlesi A.
      • Perugini J.
      • Manieri M.
      • Ciarmela P.
      • Mondini E.
      • et al.
      Molecular aspects of adipoepithelial transdifferentiation in mouse mammary gland.
      ,
      • Hennighausen L.
      • Robinson G.W.
      Information networks in the mammary gland.
      ). Prior studies have demonstrated that preadipocytes labeled with PDGFRA transdifferentiated into MECs (
      • Joshi P.A.
      • Waterhouse P.D.
      • Kasaian K.
      • Fang H.
      • Gulyaeva O.
      • Sul H.S.
      • et al.
      PDGFRalpha(+) stromal adipocyte progenitors transition into epithelial cells during lobulo-alveologenesis in the murine mammary gland.
      ). In mammary gland, the interactions between cell–cell and cell–ECM provide a complex biochemical signal network, which is very important for the process of transdifferentiation (
      • Gouon-Evans V.
      • Pollard J.W.
      Unexpected deposition of brown fat in mammary gland during postnatal development.
      ,
      • Tong J.
      • Mou S.
      • Xiong L.
      • Wang Z.
      • Wang R.
      • Weigand A.
      • et al.
      Adipose-derived mesenchymal stem cells formed acinar-like structure when stimulated with breast epithelial cells in three-dimensional culture.
      ). To study the influence of human breast cell HBL-100 on adipose-derived stem cells (ASCs) in Matrigel, the expanded ASC tended to contract and shrank lasted for 72 h when treated with the CM made from HBL-100 culture medium; and ASCs formed acinar structures with high level of KRT18 in Matrigel (
      • Tong J.
      • Mou S.
      • Xiong L.
      • Wang Z.
      • Wang R.
      • Weigand A.
      • et al.
      Adipose-derived mesenchymal stem cells formed acinar-like structure when stimulated with breast epithelial cells in three-dimensional culture.
      ). Similar to earlier findings, in our study, GM-preadipocytes expressed PDGFRA, PDGFRA (+) preadipocytes were induced by CM to transdifferentiate into MECs with high KRT18 expression, and subcellular structures with MEC features appeared. During transdifferentiation process of GM-preadipocytes, the GM-preadipocytes were cultured on Matrigel in CM for 14 days. The antennae of transdifferentiated GM-preadipocytes shrank at day 4, and the transdifferentiated GM-preadipocytes formed acinar-like structures at day 8, and there was no change from day 8 to day 14. The ultrastructural examination highlighted the existence of the adipose plasticity in adipocyte to MEC transdifferentiation (
      • De Matteis R.
      • Zingaretti M.C.
      • Murano I.
      • Vitali A.
      • Frontini A.
      • Giannulis I.
      • et al.
      In vivo physiological transdifferentiation of adult adipose cells.
      ). In line with the previous studies, our results achieved another adipose plasticity, which induced the GM-preadipocyte transdifferentiation into MECs.
      Whether the process is reversible has not been further explored. For process opposite to the transdifferentiation of preadipocytes into MECs, a prior study has proved that mouse MECs could reversibly transdifferentiate into adipocytes during mammary gland involution (
      • Giordano A.
      • Perugini J.
      • Kristensen D.M.
      • Sartini L.
      • Frontini A.
      • Kajimura S.
      • et al.
      Mammary alveolar epithelial cells convert to brown adipocytes in post-lactating mice.
      ); however, the molecular mechanism of MEC transdifferentiation in vitro has not been reported. The objectives of this work were to achieve the transdifferentiation of GM-preadipocytes into MECs with acinar features induced by the CM in vitro and identified potential candidate genes and regulatory networks that are related to the transdifferentiation, explaining why mammary adipose tissues disappear during pregnancy and lactation (
      • Cinti S.
      Pink adipocytes.
      ).
      The cell–cell interactions and paracrine connections between MECs and stroma cells in mammary gland may reveal a possible mechanism of transdifferentiation (
      • Howlett A.R.
      • Bissell M.J.
      The influence of tissue microenvironment (stroma and extracellular matrix) on the development and function of mammary epithelium.
      ,
      • McCave E.J.
      • Cass C.A.
      • Burg K.J.
      • Booth B.W.
      The normal microenvironment directs mammary gland development.
      ). Currently, the multiomics of GM-preadipocyte-to-MEC transdifferentiation was rarely examined, but we also compared the transcriptome and DNA methylome profiles between GM-preadipocytes and transdifferentiated GM-preadipocytes. The candidate genes related to the different regulatory processes of transdifferentiation were screened from DEGs and DMGs. Among these genes, DEGs and DMGs were mainly associated with immune response, cell adhesion, cell cycle, morphogenesis, and organization. In addition, phospholipase D signaling pathway and lysine degradation pathway were observed in enriched pathways of DMGs. Phospholipase D signaling pathway has participation in maintenance of cell structure and stability (
      • Gomez-Cambronero J.
      The exquisite regulation of PLD2 by a wealth of interacting proteins: S6K, Grb2, Sos, WASp and Rac2 (and a surprise discovery: PLD2 is a GEF).
      ). While lysine degradation pathway includes two aspects: (i) lysine decomposition into the citric acid and (ii) carnitine. Carnitine has been shown to be involved in milk protein synthesis in previous studies (
      • Ramanau A.
      • Kluge H.
      • Spilke J.
      • Eder K.
      Supplementation of sows with L-carnitine during pregnancy and lactation improves growth of the piglets during the suckling period through increased milk production.
      ,
      • Pirestani A.
      • Aghakhani M.
      The effects of rumen-protected choline and l-carnitine supplementation in the transition period on reproduction, production, and some metabolic diseases of dairy cattle.
      ). On the other hand, detected TFs were enriched in promoters of DEGs and DMGs, including CTCF, BACH2, STAT1, Sp5 TF (SP5), and nuclear factor I B (NFIB). After cells were stimulated with the exogenous additive, they might respond to stimulation initially by interferon produced by STAT1-mediated immune response (
      • Haricharan S.
      • Li Y.
      STAT signaling in mammary gland differentiation, cell survival and tumorigenesis.
      ,
      • Pausch H.
      • Emmerling R.
      • Schwarzenbacher H.
      • Fries R.
      A multi-trait meta-analysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle.
      ). Whereas, BACH2 is a regulator of fatty acid profiles in mammary gland (
      • Pegolo S.
      • Dadousis C.
      • Mach N.
      • Ramayo-Caldas Y.
      • Mele M.
      • Conte G.
      • et al.
      SNP co-association and network analyses identify E2F3, KDM5A and BACH2 as key regulators of the bovine milk fatty acid profile.
      ), CTCF regulates the localization of adhesin protein (
      • Zhao H.
      • Zhang S.
      • Wu X.
      • Pan C.
      • Li X.
      • Lei C.
      • et al.
      DNA methylation pattern of the goat PITX1 gene and its effects on milk performance.
      ), and NFIB regulates the formation of alveolar epithelium via the Hippo signaling pathway (
      • Grunder A.
      • Ebel T.T.
      • Mallo M.
      • Schwarzkopf G.
      • Shimizu T.
      • Sippel A.E.
      • et al.
      Nuclear factor I-B (Nfib) deficient mice have severe lung hypoplasia.
      ,
      • Gokey J.J.
      • Snowball J.
      • Sridharan A.
      • Sudha P.
      • Kitzmiller J.A.
      • Xu Y.
      • et al.
      YAP regulates alveolar epithelial cell differentiation and AGER via NFIB/KLF5/NKX2-1.
      ). These findings stated that the regulation of TFs is essential for GM-preadipocyte transdifferentiation.
      The correlation between gene expression and DNA methylation levels was also calculated. The NGs were annotated into polarity establishment, regulation of MEC adhesion and differentiation, mammary gland epithelium development, and mammary gland alveolus development. Subsequently, three potential genes that epigenetically regulate transdifferentiation were screened when overlap of DEGs and DMGs was combined with the NGs. The SPPL3 alters the pattern of cellular N-glycosylation by inducing the proteolytic release of glycosidase and glycosyltransferase extracellular domains (
      • Voss M.
      • Kunzel U.
      • Higel F.
      • Kuhn P.H.
      • Colombo A.
      • Fukumori A.
      • et al.
      Shedding of glycan-modifying enzymes by signal peptide peptidase-like 3 (SPPL3) regulates cellular N-glycosylation.
      ), and the rate of glycan flux/the level of N-glycan biosynthetic enzyme determines the cellular behavior in reprogramming (
      • Kasper D.M.
      • Hintzen J.
      • Wu Y.
      • Ghersi J.J.
      • Mandl H.K.
      • Salinas K.E.
      • et al.
      The N-glycome regulates the endothelial-to-hematopoietic transition.
      ). Silencing RND1 reduced MEC adhesion and polarity, whereas overexpression improved E-cadherin expression (
      • Qin C.D.
      • Ma D.N.
      • Zhang S.Z.
      • Zhang N.
      • Ren Z.G.
      • Zhu X.D.
      • et al.
      The Rho GTPase Rnd1 inhibits epithelial-mesenchymal transition in hepatocellular carcinoma and is a favorable anti-metastasis target.
      ). The PLCL1 could activate the browning of adipocyte plasticity (
      • Xiong Z.
      • Xiao W.
      • Bao L.
      • Xiong W.
      • Xiao H.
      • Qu Y.
      • et al.
      Tumor cell "slimming" regulates tumor progression through PLCL1/UCP1-mediated lipid browning.
      ). These reports illustrated that these candidate genes are crucial for GM-preadipocyte-to-MEC transdifferentiation. In addition, the secretory proteome in the CM was detected to construct complete regulatory networks. The epigenetic regulatory network was constructed by interacting DMGs and secretory proteins. In the regulatory network, ALB, ACTG1, NID1, and HSP90AA1 were deemed to be core proteins. Among these proteins, ALB involves in negative regulation of apoptotic process (
      • Fu X.
      • Yang Y.
      • Zhang D.
      Molecular mechanism of albumin in suppressing invasion and metastasis of hepatocellular carcinoma.
      ); ACTG1 regulates the cell adhesion, morphogenesis of a polarized epithelium, and the level of E-cadherin (
      • Huang S.
      • Cai M.
      • Zheng Y.
      • Zhou L.
      • Wang Q.
      • Chen L.
      miR-888 in MCF-7 side population sphere cells directly targets E-cadherin.
      ); NID1 promotes cell adhesion by connecting collagen IV and laminin (
      • Jagroop R.
      • Martin C.J.
      • Moorehead R.A.
      Nidogen 1 regulates proliferation and migration/invasion in murine claudin-low mammary tumor cells.
      ); and HSP90AA1 participates in cell proliferation and apoptosis (
      • Liu C.
      • Lin C.
      • Wang D.
      • Wang J.
      • Tao Y.
      • Li Y.
      • et al.
      Procr functions as a signaling receptor and is essential for the maintenance and self-renewal of mammary stem cells.
      ). In epigenetic regulatory axis, RND1 was considered as a critical candidate gene regulated by secretory protein, which needs to be further verified in subsequent experiments.
      Nevertheless, there are other limitations: the proteomes of GM-preadipocytes and transdifferentiated GM-preadipocytes were not detected, and the omics research is insufficient. Collectively in this study, three candidate genes and four candidate proteins were identified through integrated analysis of transcriptome, DNA methylome, and secretory proteome. These candidate functional factors and paracrine mediators play key regulatory roles in the whole processes of GM-preadipocyte-to-MEC transdifferentiation. This work provides new insights into the molecular mechanisms of adipocyte plasticity and mammary development during pregnancy and lactation.

      Experimental procedures

      Experimental design

      The GM-preadipocytes were isolated from adipose tissues around the goat mammary gland, whereas the GMECs were isolated from goat mammary gland. These cells were cultured in Matrigel (Corning), and the GMEC culture medium was collected as CM to induce the GM-preadipocyte transdifferentiation. Subsequently, GM-preadipocytes and transdifferentiated GM-preadipocytes were used for transcriptome sequencing and WGBS. Moreover, secretory proteins in CM were detected. The candidate genes, proteins, and epigenetic regulatory axes were screened through integration of transcriptome, DNA methylome, and secretory proteome (Fig. 8).

      Isolation, culture, and induced differentiation of GM-preadipocytes

      The animal research ethics committee of Northwest A&F University approved the animal use protocol. GM-preadipocytes were isolated from adipose tissues around the goat mammary gland under sterile conditions. The adipose tissues were cut into 1 mm3 and digested at 37 °C with 2% (w/v) collagenase type I (Invitrogen) and then diluted with Dulbecco's modified Eagle's medium (DMEM)/F12 (Hyclone) for 2 to 3 h, and the cell suspension was centrifuged at 700g for 10 min. Then, pellets were resuspended in the growth medium containing 89% DMEM/F12, 10% fetal bovine serum (FBS; Gemini) and 1% penicillin–streptomycin (pen–strep; Solarbio). The undigested parts were filtered using 40-μm cell strainer, and harvested GM-preadipocytes were cultured in growth medium at 37 °C with 5% CO2.
      Once cell confluence reached 90 to 100%, the growth-arrested cells were initiated to induce differentiation by treatment of the adipogenic induction medium containing 89% DMEM/F12, 10% FBS, 1% pen–strep, 1 μmol/l of dexamethasone (Solarbio), 250 μmol/l of 3-isobutyl-1-methylxanthine (Solarbio), and 5 μg/ml of insulin (Solarbio). After 2 days of adipogenic induction, cells were cultured in maintenance medium (89% DMEM/F12, 10% FBS, 1% pen–strep, and 5 μg/ml of insulin) for 8 days to maintain adipogenic induction, and maintenance medium was changed every 2 days.

      Isolation and culture of GMECs

      GMECs were isolated from healthy peak-lactation goats as described previously (
      • Shi H.
      • Shi H.
      • Luo J.
      • Wang W.
      • Haile A.B.
      • Xu H.
      • et al.
      Establishment and characterization of a dairy goat mammary epithelial cell line with human telomerase (hT-MECs).
      ). Briefly, mammary tissue was collected and rinsed repeatedly with D-Hank's (Hyclone) containing 3% pen–strep to remove adipose tissue and connective tissue. Parenchyma tissue was cut into 1 mm3 and placed in petri dishes pretreated with FBS and a few drops of DMEM/F12 growth medium consisting of 10% FBS, 1% pen–strep, 1% insulin, 0.1% epidermal growth factor (Invitrogen), and 2% hydrocortisone (Sigma). Near the tissue block, added a few drops of growth medium every 30 min to keep the tissue block wet. Two hours later, 3 ml of growth medium was added for tissue block culture, and GMECs were isolated and purified by differential adhesion. When cells were grown to 80 to 90% confluence, the medium was discarded and cells were treated with 0.25% trypsin–EDTA solution (Solarbio) for 5 min. Afterward, the collected cell suspension was centrifuged at 1000g for 4 min, and harvested cells were cultured in the growth medium at 37 °C with 5% CO2. The growth medium was changed every 2 days.

      Oil Red O staining

      The Oil Red O staining was used to detect changes in lipid droplets in GM-preadipocytes at day 2, 5, and 8 of adipogenic induction. Briefly, cells were washed three times with PBS (Hyclone) and fixed with 4% paraformaldehyde (Solarbio) overnight at room temperature. Then, cells were dyed with working solution of Oil Red O (Solarbio), and dye was extracted by 60% isopropanol for 5 min at room temperature. The images were captured using an inverted microscope.

      Immunofluorescence

      To verify some of the markers for GM-preadipocytes and GMECs at the protein level, we analyzed GM-preadipocytes markers (DLK1 and CD34) and GMEC marker (KRT18) by antibody-based immunofluorescence staining. The steps are as follows: at 90 to 100% cell confluence, then cells were washed three times with PBS containing 1% Tween-20 (Solarbio) and fixed with 4% paraformaldehyde. After 30 min, cells were treated with 0.1% Triton X-100 (Solarbio) and 5% bovine serum albumin (Solarbio) for 10 min, respectively. Afterward, cells were incubated with primary antibodies (Table S17) at 4 °C overnight on a low-speed shaker and then washed in PBS containing 1% Tween-20. After this step, FITC-conjugated secondary antibodies (Proteintech) were added, and cells were incubated for 1 h at room temperature. The nucleus was counterstained with 4,6-diamino-2-phenyl indole (Solarbio) for 5 min. Finally, the cells were visualized under the fluorescence microscope and confocal laser scanning microscope.

      Matrigel coating and cell culturing

      The GM-preadipocytes and GMECs were cultured in Matrigel as described previously (
      • Muthuswamy S.K.
      • Li D.
      • Lelievre S.
      • Bissell M.J.
      • Brugge J.S.
      ErbB2, but not ErbB1, reinitiates proliferation and induces luminal repopulation in epithelial acini.
      ,
      • Debnath J.
      • Muthuswamy S.K.
      • Brugge J.S.
      Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures.
      ). Matrigel was thawed in ice and diluted to 3 mg/ml with DMEM/F12. The bottoms of 48-well plate were coated with 250 μl Matrigel, and Matrigel-coated plates were placed at 37 °C for 1 h to solidify. Two types of cells were suspended using growth medium with 2% Matrigel (3D growth medium), inoculated into Matrigel-coated plates, and cultured in incubator at 37 °C with 5% CO2, respectively. The morphological changes were observed at the live cell imaging system. The medium was carefully changed every 2 days.

      Preparation of CM

      During the process of CM preparation, GMECs were continuously cultured in Matrigel. The GMECs were cultured in Matrigel for 2 days and washed three times with preheated PBS. In the following 2-day culture, the 3D growth medium was replaced by serum-free 3D growth medium for culture medium collection. The resulting medium was centrifuged at 300g for 5 min, and then supernatant was collected and filtered with 0.22-μm filter. The filtered supernatant mixed with an equal volume of DMEM/F12 was subjected to a quality inspection, which included sterility testing and appearance inspection (
      • Tong J.
      • Mou S.
      • Xiong L.
      • Wang Z.
      • Wang R.
      • Weigand A.
      • et al.
      Adipose-derived mesenchymal stem cells formed acinar-like structure when stimulated with breast epithelial cells in three-dimensional culture.
      ,
      • Pallegar N.K.
      • Garland C.J.
      • Mahendralingam M.
      • Viloria-Petit A.M.
      • Christian S.L.
      A novel 3-dimensional co-culture method reveals a partial mesenchymal to epithelial transition in breast cancer cells induced by adipocytes.
      ). Finally, qualified medium was used as CM in induction of transdifferentiation and protein mass spectrometry.

      Induction of transdifferentiation

      The transdifferentiation of GM-preadipocyte into MEC was carried out in Matrigel. The GM-preadipocytes were cultured for 24 h with serum-free 3D growth medium. After 24 h, cells were washed three times with preheated PBS, and the preheated CM was used to induce the GM-preadipocyte transdifferentiation (Gpreadipocyte group). Simultaneously, GM-preadipocytes were cultured with the serum-free 3D growth medium in Matrigel as the control group (Preadipocyte group). The medium was changed every 2 days, and induction of transdifferentiation was kept for 14 days. The morphological changes were observed every day using live cell imaging system.

      TEM

      TEM was used to observe subcellular structures of GM-preadipocytes and transdifferentiated GM-preadipocytes as described previously (
      • Ge T.
      • Ye Y.
      • Zhang H.
      Ultrastructure of telocytes, a new type of interstitial cells in the myocardium of the Chinese giant salamander (Andrias davidianus).
      ). For cells of Preadipocyte and Gpreadipocyte groups, the medium was removed at day 8, and cells were washed three times with preheated PBS. Subsequently, the cells were fixed with 2.5% glutaraldehyde (0.1 M, pH = 7.2–7.4) for 5 h and washed three times with 0.1 M phosphoric acid buffer. Then, cells were fixed with 1% osmic acid (0.1 M, pH = 7.2–7.4) for 2 to 3 h and washed three times with 0.1 M phosphoric acid buffer for 15 min. Followed by twice gradient dehydration in a concentration series of ethanol (30, 50, 70, 80, 90, and 100% ethanol for 10 min, respectively). The samples were infiltrated with resin and embedded in paraffin using a capsule tool. Then, samples were fixed in an incubator at 55 °C for 48 h and cut into 50 to 70 nm-thick slices using an ultrathin slicer. Finally, samples were stained with 2% uranyl acetate and lead citrate for 30 and 15 min, respectively, and observed under TEM.

      Transcriptome data

      Sample collection and RNA preparation

      GM-preadipocytes and transdifferentiated GM-preadipocytes were collected to perform transcriptome sequencing. Cell recovery solution (Corning) was used to recover cells from Matrigel according to the manufacturer's instructions. The cell culture dishes were incubated with prechilled cell recovery solution at 4 °C for 20 min. After this step, the cells were separated from cell recovery solution by brief centrifugation (700g for 10 min). Finally, total RNA was extracted from collected cells using Trizol reagent (TaKaRa) according to the manufacturer’s protocol. The integrity of RNA was further detected by Agilent 2100 bioanalyzer (Agilent Technologies).

      Complementary DNA library generation, sequencing, and mapping

      The NEBNext Ultra RNA Library Prep Kit for Illumina (
      • Parkhomchuk D.
      • Borodina T.
      • Amstislavskiy V.
      • Banaru M.
      • Hallen L.
      • Krobitsch S.
      • et al.
      Transcriptome analysis by strand-specific sequencing of complementary DNA.
      ) was used for complementary DNA (cDNA) library generation. The mRNA with ployA was enriched by oligo(dT) magnetic beads and further fragmented in NEB fragmentation buffer. First-strand cDNA was synthesized in M-MuLV reverse transcriptase system, and second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The 250 to 300 bp cDNA was screened for PCR amplification, PCR product was purified using AMPure XP beads, and preliminary quality of library was assessed on the Qubit 2.0 Fluorometer. Finally, the library was diluted to 1.5 ng/μl, and insert size of library was detected by Agilent 2100 bioanalyzer. The final libraries were subjected to sequence using Illumina NovaSeq 6000 (Illumina). The obtained raw reads were cleaned to remove adapter and low-quality reads using fastp (
      • Chen S.
      • Zhou Y.
      • Chen Y.
      • Gu J.
      fastp: an ultra-fast all-in-one FASTQ preprocessor.
      ) (version 0.19.7; the parameters: -q 5 -u 50 -n 15 -l 150). Then, HISAT2 (https://github.com/DaehwanKimLab/hisat2) (
      • Kim D.
      • Langmead B.
      • Salzberg S.L.
      HISAT: a fast spliced aligner with low memory requirements.
      ) was used to map clean reads to goat reference genome ARS1. HISAT2 generated alignment was subjected to Subread featureCounts tool for generating a matrix of genewise counts.

      Identification and functional enrichment analysis of DEGs

      The DEGs were analyzed using the R package “DESeq2” (
      • Love M.I.
      • Huber W.
      • Anders S.
      Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
      ), which is based on read counts for difference analysis. |Log2(foldchange)| ≥ 1 and p < 0.05 were used as the screening criteria. R package “piano” (
      • Varemo L.
      • Nielsen J.
      • Nookaew I.
      Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods.
      ) was used to distinguish the function of downregulated and upregulated DEGs by gene set enrichment analysis, and the consensusScores function in “Piano” was used to estimate consensus gene set scores for each directionality class based on findings (gene set p values).

      RT–qPCR assay and statistical analysis

      Total RNA was reverse-transcribed into cDNA using PrimerScript RT Reagent Kit (Takara). The GAPDH was used as an internal control for gene expression analysis, and the primers used for RT–qPCR are listed in Table S18. The RT–qPCR was performed on the CFX Manager 3.1 analyzer system (Bio-Rad) using SYBR kit (Takara), and process was run as follows: 95 °C for 5 min, followed by 39 cycles of 95 °C for 30 s, 63 °C for 30 s, and 72 °C for 1 min. The expressions of genes were calculated using 2−△△Cq method. The data of RT–qPCR were analyzed using unpaired t test in GraphPad Prism 6.0 software (GraphPad Software, Inc). For all analyses, p < 0.05 was considered statistically significant.

      WGBS data

      Sample collection and DNA preparation

      GM-preadipocytes and transdifferentiated GM-preadipocytes were collected for WGBS. The total DNA was extracted using HiPure Tissue DNA Micro Kit (Magnetec) according to the manufacturer’s instructions.

      Library preparation, sequencing, mapping, and methylation calling

      The DNA quality was detected using Qubit (Thermo), and the WGBS library was prepared as described previously (
      • Zhang S.
      • Qin C.
      • Cao G.
      • Guo L.
      • Feng C.
      • Zhang W.
      Genome-wide analysis of DNA methylation profiles in a senescence-accelerated mouse prone 8 brain using whole-genome bisulfite sequencing.
      ). Covaris E220 focused-ultra sonicator was used for DNA fragmentation to produce DNA fragments ranging from 300 to 700 bp. Followed by a terminal repair and adenylation reaction, the cytosine methylation barcode was attached to DNA fragment. These DNA fragments were treated twice with bisulfite using EZ DNA Methylation-Gold Kit (Zymo Research). The Agilent Bioanalyzer 2100 system was then used to measure the insert size, and sequencing was carried out using DNBSEQ Platform (MGI Tech Co, Ltd). Finally, 100 bp paired-end reads were generated. The low-quality reads in raw data were filtered using fastp with parameter: -q 40 -u 20 -n 5 -l 30 -p 20 -w 4. BiSulfite Bolt (https://github.com/NuttyLogic/BSBolt) (
      • Farrell C.
      • Thompson M.
      • Tosevska A.
      • Oyetunde A.
      • Pellegrini M.
      BiSulfite Bolt: a bisulfite sequencing analysis platform.
      ) was carried out to map clean reads to goat reference genome ARS1, and matrices containing methylation values and counts of methylated and total bases at each site were outputted for following analysis.

      Identification of DMLs and DMRs

      Based on methylation site information, the R package “Methylkit” (
      • Akalin A.
      • Kormaksson M.
      • Li S.
      • Garrett-Bakelman F.E.
      • Figueroa M.E.
      • Melnick A.
      • et al.
      methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles.
      ) was used to analyze the DMLs and DMRs with 400-bp sliding window and 400-bp step. The DMR was defined as a region containing ≥4 DMLs. In addition, Methylkit was also used to annotate DMLs and DMRs obtaining DMGs.

      TF-binding motif enrichment analysis

      TF-binding motif enrichment analysis of promoter of DEGs and DMGs was carried out by MEME of JASPAR (
      • Amith M.T.
      • Fujimoto K.
      • Tao C.
      NET-EXPO: a gephi plugin towards social network analysis of network exposure for unipartite and bipartite graphs.
      ,
      • Fornes O.
      • Castro-Mondragon J.A.
      • Khan A.
      • van der Lee R.
      • Zhang X.
      • Richmond P.A.
      • et al.
      Jaspar 2020: update of the open-access database of transcription factor binding profiles.
      ). Promoters were defined as 5000-bp upstream and 5000-bp downstream from the transcriptional start site of each gene.

      Proteome data

      Sample preparation

      The collected CM for transdifferentiation induction was centrifuged with 200g for 10 min at 4 °C and then filtered with a 0.22-μm filter and added protease inhibitors (Roche). The resulting solution was dried by vacuum centrifugation, and it was redissolved with 1 ml of 50 mM NH4HCO3. For ultrafiltration, the Amicon Ultra-0.5 centrifugal filter device (Millipore) was pretreated according to the manufacturer's instructions. Then, the solution was concentrated using Amicon Ultra-0.5 centrifugal filters with 3 kDa cutoff (Merck Millipore). Filter-aided sample preparation was performed as described previously (
      • Wisniewski J.R.
      • Zougman A.
      • Nagaraj N.
      • Mann M.
      Universal sample preparation method for proteome analysis.
      ). Briefly, DTTs were added to the protein sample as reducing agent and then incubated at 37 °C for 1 h, and iodoacetamide was added as alkylating agent and reacted for 30 min at room temperature in the dark. After incubation, 0.5 μg of trypsin was added for an overnight digestion at 37 °C with shaking. The peptides were finally eluted and desalted on 100 μl C18 column (Thermo). Then, collected peptides were dried by vacuum centrifugation and redissolved with 10 μl 0.1% methanoic acid.

      LC–MS/MS acquisition

      The final solution was subjected to LC–MS/MS analysis using Orbitrap fusion mass spectrometer (Thermo). For LC–MS/MS analysis, using 0.1% methanoic acid and 0.1% methanoic acid, and 80% acetonitrile as mobile phase A and mobile phase B, respectively; The flow rate was adjusted to 300 nl/min; and analysis time, 0 min (8% phase B), 38 min (38% phase B), 43 min (48% phase B), 45 min (100% phase B), and 55 min (100% phase B). Default parameters were used for MS scan. We initially used Unified Protein Database (UniProt; https://www.uniprot.org/) to create a sample-specific database for identifying goat protein sequences. The Orbi RAW files were searched directly using Byonic (Protein Metrics) on Proteome Discoverer, version 2.2 software (Thermo), and spectra were searched by sample-specific database. The cleavage site of trypsin, the ion mass tolerance of 0.8 Da fragment, and the parent ion tolerance of 10.0 ppm were used as parameters. In addition, cysteine carbamidomethylation was defined as fixed and methionine oxidation and N-terminal acetylation as variable modifications. The proteome with high confidence and at least one unique peptide was collected for further analysis.

      Functional enrichment analysis

      DAVID (https://david.ncifcrf.gov/) was used to perform GO enrichment analysis, and KOBAS (http://kobas.cbi.pku.edu.cn/kobas3) was used for KEGG pathway functional enrichment analysis. These results were visualized using R package “ggplot2”.

      Data availability

      All data are contained within the article, except the data deposited into a publicly accessible repository as follows:
      The data were deposited in China National GeneBank (https://db.cngb.org/search/project/CNP0003568/ and https://db.cngb.org/search/project/CNP0003010/) embargoed until publication.

      Supporting information

      This article contains supporting information.

      Conflict of interest

      The authors declare that they have no conflicts of interest with the contents of this article.

      Acknowledgments

      We thank the support of China National GeneBank. This work was supported by the National Natural Science Foundation of China (grant no.: 31572368 ) and the Key Projects of Shaanxi Natural Science Foundation Research Plan (grant no.: 2022JZ-10 ).

      Author contributions

      X.-R. Y. and H.-L. Z. conceptualization; X.-R. Y. and T. S. software; H.-L. Z. validation; X.-R. Y., T. S., and Y.-F. L. formal analysis; X.-R. Y., T. S., J.-Y. X., and Y.-F. L. investigation; X.-R. Y. and T. S. data curation; X.-R. Y. writing–original draft; X.-R. Y. and H.-L. Z. writing–review & editing; X.-R. Y. and T. S. visualization; H.-L. Z. supervision; H.-L. Z. project administration; H.-L. Z. funding acquisition.

      References

        • Neville M.C.
        • Medina D.
        • Monks J.
        • Hovey R.C.
        The mammary fat pad.
        J. Mammary Gland Biol. Neoplasia. 1998; 3: 109-116
        • Giordano A.
        • Perugini J.
        • Kristensen D.M.
        • Sartini L.
        • Frontini A.
        • Kajimura S.
        • et al.
        Mammary alveolar epithelial cells convert to brown adipocytes in post-lactating mice.
        J. Cell Physiol. 2017; 232: 2923-2928
        • Giordano A.
        • Smorlesi A.
        • Frontini A.
        • Barbatelli G.
        • Cinti S.
        White, brown and pink adipocytes: the extraordinary plasticity of the adipose organ.
        Eur. J. Endocrinol. 2014; 170: R159-171
        • Smith G.H.
        • Medina D.
        Re-evaluation of mammary stem cell biology based on in vivo transplantation.
        Breast Cancer Res. 2008; 10: 203
        • Joshi P.A.
        • Waterhouse P.D.
        • Kasaian K.
        • Fang H.
        • Gulyaeva O.
        • Sul H.S.
        • et al.
        PDGFRalpha(+) stromal adipocyte progenitors transition into epithelial cells during lobulo-alveologenesis in the murine mammary gland.
        Nat. Commun. 2019; 10: 1760
        • Colleluori G.
        • Perugini J.
        • Barbatelli G.
        • Cinti S.
        Mammary gland adipocytes in lactation cycle, obesity and breast cancer.
        Rev. Endocr. Metab. Disord. 2021; 22: 241-255
        • Wang Q.A.
        • Song A.
        • Chen W.
        • Schwalie P.C.
        • Zhang F.
        • Vishvanath L.
        • et al.
        Reversible de-differentiation of mature white adipocytes into preadipocyte-like precursors during lactation.
        Cell Metab. 2018; 28: 282-288.e283
        • Lee Y.H.
        • Petkova A.P.
        • Mottillo E.P.
        • Granneman J.G.
        In vivo identification of bipotential adipocyte progenitors recruited by beta3-adrenoceptor activation and high-fat feeding.
        Cell Metab. 2012; 15: 480-491
        • Carriere A.
        • Jeanson Y.
        • Cote J.A.
        • Dromard C.
        • Galinier A.
        • Menzel S.
        • et al.
        Identification of the ectoenzyme CD38 as a marker of committed preadipocytes.
        Int. J. Obes. (Lond). 2017; 41: 1539-1546
        • Berry R.
        • Rodeheffer M.S.
        Characterization of the adipocyte cellular lineage in vivo.
        Nat. Cell Biol. 2013; 15: 302-308
        • Faridi N.
        • Bathaie S.Z.
        • Abroun S.
        • Farzaneh P.
        • Karbasian H.
        • Tamanoi F.
        • et al.
        Isolation and characterization of the primary epithelial breast cancer cells and the adjacent normal epithelial cells from Iranian women's breast cancer tumors.
        Cytotechnology. 2018; 70: 625-639
        • Weng Y.R.
        • Cui Y.
        • Fang J.Y.
        Biological functions of cytokeratin 18 in cancer.
        Mol. Cancer Res. 2012; 10: 485-493
        • Hovey R.C.
        • Trott J.F.
        Morphogenesis of mammary gland development.
        Adv. Exp. Med. Biol. 2004; 554: 219-228
        • Yin Y.
        • Yuan H.
        • Wang C.
        • Pattabiraman N.
        • Rao M.
        • Pestell R.G.
        • et al.
        3-phosphoinositide-dependent protein kinase-1 activates the peroxisome proliferator-activated receptor-gamma and promotes adipocyte differentiation.
        Mol. Endocrinol. 2006; 20: 268-278
        • Whyte J.
        • Thornton L.
        • McNally S.
        • McCarthy S.
        • Lanigan F.
        • Gallagher W.M.
        • et al.
        PKCzeta regulates cell polarisation and proliferation restriction during mammary acinus formation.
        J. Cell Sci. 2010; 123: 3316-3328
        • Muthuswamy S.K.
        • Li D.
        • Lelievre S.
        • Bissell M.J.
        • Brugge J.S.
        ErbB2, but not ErbB1, reinitiates proliferation and induces luminal repopulation in epithelial acini.
        Nat. Cell Biol. 2001; 3: 785-792
        • Gusterson B.A.
        • Stein T.
        Human breast development.
        Semin. Cell Dev. Biol. 2012; 23: 567-573
        • Nagy D.
        • Gillis C.M.C.
        • Davies K.
        • Fowden A.L.
        • Rees P.
        • Wills J.W.
        • et al.
        Developing ovine mammary terminal duct lobular units have a dynamic mucosal and stromal immune microenvironment.
        Commun. Biol. 2021; 4: 993
        • Lee E.Y.
        • Lee W.H.
        • Kaetzel C.S.
        • Parry G.
        • Bissell M.J.
        Interaction of mouse mammary epithelial cells with collagen substrata: regulation of casein gene expression and secretion.
        Proc. Natl. Acad. Sci. U. S. A. 1985; 82: 1419-1423
        • Sumbal J.
        • Chiche A.
        • Charifou E.
        • Koledova Z.
        • Li H.
        Primary mammary organoid model of lactation and involution.
        Front. Cell Dev. Biol. 2020; 8: 68
        • Mackenzie D.D.
        • Forsyth I.A.
        • Brooker B.E.
        • Turvey A.
        Culture of bovine mammary epithelial cells on collagen gels.
        Tissue Cell. 1982; 14: 231-241
        • Riley L.G.
        • Gardiner-Garden M.
        • Thomson P.C.
        • Wynn P.C.
        • Williamson P.
        • Raadsma H.W.
        • et al.
        The influence of extracellular matrix and prolactin on global gene expression profiles of primary bovine mammary epithelial cells in vitro.
        Anim. Genet. 2010; 41: 55-63
        • Kozlowski M.
        • Gajewska M.
        • Majewska A.
        • Jank M.
        • Motyl T.
        Differences in growth and transcriptomic profile of bovine mammary epithelial monolayer and three-dimensional cell cultures.
        J. Physiol. Pharmacol. 2009; 60: 5-14
        • Wang H.
        • Lacoche S.
        • Huang L.
        • Xue B.
        • Muthuswamy S.K.
        Rotational motion during three-dimensional morphogenesis of mammary epithelial acini relates to laminin matrix assembly.
        Proc. Natl. Acad. Sci. U. S. A. 2013; 110: 163-168
        • Cinti S.
        Transdifferentiation properties of adipocytes in the adipose organ.
        Am. J. Physiol. Endocrinol. Metab. 2009; 297: E977-986
        • Seale P.
        • Conroe H.M.
        • Estall J.
        • Kajimura S.
        • Frontini A.
        • Ishibashi J.
        • et al.
        Prdm16 determines the thermogenic program of subcutaneous white adipose tissue in mice.
        J. Clin. Invest. 2011; 121: 96-105
        • Morroni M.
        • Giordano A.
        • Zingaretti M.C.
        • Boiani R.
        • De Matteis R.
        • Kahn B.B.
        • et al.
        Reversible transdifferentiation of secretory epithelial cells into adipocytes in the mammary gland.
        Proc. Natl. Acad. Sci. U. S. A. 2004; 101: 16801-16806
        • De Matteis R.
        • Zingaretti M.C.
        • Murano I.
        • Vitali A.
        • Frontini A.
        • Giannulis I.
        • et al.
        In vivo physiological transdifferentiation of adult adipose cells.
        Stem Cells. 2009; 27: 2761-2768
        • Prokesch A.
        • Smorlesi A.
        • Perugini J.
        • Manieri M.
        • Ciarmela P.
        • Mondini E.
        • et al.
        Molecular aspects of adipoepithelial transdifferentiation in mouse mammary gland.
        Stem Cells. 2014; 32: 2756-2766
        • Hennighausen L.
        • Robinson G.W.
        Information networks in the mammary gland.
        Nat. Rev. Mol. Cell Biol. 2005; 6: 715-725
        • Gouon-Evans V.
        • Pollard J.W.
        Unexpected deposition of brown fat in mammary gland during postnatal development.
        Mol. Endocrinol. 2002; 16: 2618-2627
        • Tong J.
        • Mou S.
        • Xiong L.
        • Wang Z.
        • Wang R.
        • Weigand A.
        • et al.
        Adipose-derived mesenchymal stem cells formed acinar-like structure when stimulated with breast epithelial cells in three-dimensional culture.
        PLoS One. 2018; 13e0204077
        • Cinti S.
        Pink adipocytes.
        Trends Endocrinol. Metab. 2018; 29: 651-666
        • Howlett A.R.
        • Bissell M.J.
        The influence of tissue microenvironment (stroma and extracellular matrix) on the development and function of mammary epithelium.
        Epithelial Cell Biol. 1993; 2: 79-89
        • McCave E.J.
        • Cass C.A.
        • Burg K.J.
        • Booth B.W.
        The normal microenvironment directs mammary gland development.
        J. Mammary Gland Biol. Neoplasia. 2010; 15: 291-299
        • Gomez-Cambronero J.
        The exquisite regulation of PLD2 by a wealth of interacting proteins: S6K, Grb2, Sos, WASp and Rac2 (and a surprise discovery: PLD2 is a GEF).
        Cell Signal. 2011; 23: 1885-1895
        • Ramanau A.
        • Kluge H.
        • Spilke J.
        • Eder K.
        Supplementation of sows with L-carnitine during pregnancy and lactation improves growth of the piglets during the suckling period through increased milk production.
        J. Nutr. 2004; 134: 86-92
        • Pirestani A.
        • Aghakhani M.
        The effects of rumen-protected choline and l-carnitine supplementation in the transition period on reproduction, production, and some metabolic diseases of dairy cattle.
        J. Appl. Anim. Res. 2018; 46: 435-440
        • Haricharan S.
        • Li Y.
        STAT signaling in mammary gland differentiation, cell survival and tumorigenesis.
        Mol. Cell Endocrinol. 2014; 382: 560-569
        • Pausch H.
        • Emmerling R.
        • Schwarzenbacher H.
        • Fries R.
        A multi-trait meta-analysis with imputed sequence variants reveals twelve QTL for mammary gland morphology in Fleckvieh cattle.
        Genet. Sel Evol. 2016; 48: 14
        • Pegolo S.
        • Dadousis C.
        • Mach N.
        • Ramayo-Caldas Y.
        • Mele M.
        • Conte G.
        • et al.
        SNP co-association and network analyses identify E2F3, KDM5A and BACH2 as key regulators of the bovine milk fatty acid profile.
        Sci. Rep. 2017; 717317
        • Zhao H.
        • Zhang S.
        • Wu X.
        • Pan C.
        • Li X.
        • Lei C.
        • et al.
        DNA methylation pattern of the goat PITX1 gene and its effects on milk performance.
        Arch. Anim. Breed. 2019; 62: 59-68
        • Grunder A.
        • Ebel T.T.
        • Mallo M.
        • Schwarzkopf G.
        • Shimizu T.
        • Sippel A.E.
        • et al.
        Nuclear factor I-B (Nfib) deficient mice have severe lung hypoplasia.
        Mech. Dev. 2002; 112: 69-77
        • Gokey J.J.
        • Snowball J.
        • Sridharan A.
        • Sudha P.
        • Kitzmiller J.A.
        • Xu Y.
        • et al.
        YAP regulates alveolar epithelial cell differentiation and AGER via NFIB/KLF5/NKX2-1.
        iScience. 2021; 24102967
        • Voss M.
        • Kunzel U.
        • Higel F.
        • Kuhn P.H.
        • Colombo A.
        • Fukumori A.
        • et al.
        Shedding of glycan-modifying enzymes by signal peptide peptidase-like 3 (SPPL3) regulates cellular N-glycosylation.
        EMBO J. 2014; 33: 2890-2905
        • Kasper D.M.
        • Hintzen J.
        • Wu Y.
        • Ghersi J.J.
        • Mandl H.K.
        • Salinas K.E.
        • et al.
        The N-glycome regulates the endothelial-to-hematopoietic transition.
        Science. 2020; 370: 1186-1191
        • Qin C.D.
        • Ma D.N.
        • Zhang S.Z.
        • Zhang N.
        • Ren Z.G.
        • Zhu X.D.
        • et al.
        The Rho GTPase Rnd1 inhibits epithelial-mesenchymal transition in hepatocellular carcinoma and is a favorable anti-metastasis target.
        Cell Death Dis. 2018; 9: 486
        • Xiong Z.
        • Xiao W.
        • Bao L.
        • Xiong W.
        • Xiao H.
        • Qu Y.
        • et al.
        Tumor cell "slimming" regulates tumor progression through PLCL1/UCP1-mediated lipid browning.
        Adv. Sci. (Weinh). 2019; 61801862
        • Fu X.
        • Yang Y.
        • Zhang D.
        Molecular mechanism of albumin in suppressing invasion and metastasis of hepatocellular carcinoma.
        Liver Int. 2021; 42: 696-709
        • Huang S.
        • Cai M.
        • Zheng Y.
        • Zhou L.
        • Wang Q.
        • Chen L.
        miR-888 in MCF-7 side population sphere cells directly targets E-cadherin.
        J. Genet. Genomics. 2014; 41: 35-42
        • Jagroop R.
        • Martin C.J.
        • Moorehead R.A.
        Nidogen 1 regulates proliferation and migration/invasion in murine claudin-low mammary tumor cells.
        Oncol. Lett. 2021; 21: 52
        • Liu C.
        • Lin C.
        • Wang D.
        • Wang J.
        • Tao Y.
        • Li Y.
        • et al.
        Procr functions as a signaling receptor and is essential for the maintenance and self-renewal of mammary stem cells.
        Cell Rep. 2022; 38110548
        • Shi H.
        • Shi H.
        • Luo J.
        • Wang W.
        • Haile A.B.
        • Xu H.
        • et al.
        Establishment and characterization of a dairy goat mammary epithelial cell line with human telomerase (hT-MECs).
        Anim. Sci. J. 2014; 85: 735-743
        • Debnath J.
        • Muthuswamy S.K.
        • Brugge J.S.
        Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures.
        Methods. 2003; 30: 256-268
        • Pallegar N.K.
        • Garland C.J.
        • Mahendralingam M.
        • Viloria-Petit A.M.
        • Christian S.L.
        A novel 3-dimensional co-culture method reveals a partial mesenchymal to epithelial transition in breast cancer cells induced by adipocytes.
        J. Mammary Gland Biol. Neoplasia. 2019; 24: 85-97
        • Ge T.
        • Ye Y.
        • Zhang H.
        Ultrastructure of telocytes, a new type of interstitial cells in the myocardium of the Chinese giant salamander (Andrias davidianus).
        Eur. J. Histochem. 2019; 63: 3021
        • Parkhomchuk D.
        • Borodina T.
        • Amstislavskiy V.
        • Banaru M.
        • Hallen L.
        • Krobitsch S.
        • et al.
        Transcriptome analysis by strand-specific sequencing of complementary DNA.
        Nucl. Acids Res. 2009; 37: e123
        • Chen S.
        • Zhou Y.
        • Chen Y.
        • Gu J.
        fastp: an ultra-fast all-in-one FASTQ preprocessor.
        Bioinformatics. 2018; 34: i884-i890
        • Kim D.
        • Langmead B.
        • Salzberg S.L.
        HISAT: a fast spliced aligner with low memory requirements.
        Nat Methods. 2015; 12: 357-360
        • Love M.I.
        • Huber W.
        • Anders S.
        Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
        Genome Biol. 2014; 15: 550
        • Varemo L.
        • Nielsen J.
        • Nookaew I.
        Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods.
        Nucl. Acids Res. 2013; 41: 4378-4391
        • Zhang S.
        • Qin C.
        • Cao G.
        • Guo L.
        • Feng C.
        • Zhang W.
        Genome-wide analysis of DNA methylation profiles in a senescence-accelerated mouse prone 8 brain using whole-genome bisulfite sequencing.
        Bioinformatics. 2017; 33: 1591-1595
        • Farrell C.
        • Thompson M.
        • Tosevska A.
        • Oyetunde A.
        • Pellegrini M.
        BiSulfite Bolt: a bisulfite sequencing analysis platform.
        Gigascience. 2021; 10: giab033
        • Akalin A.
        • Kormaksson M.
        • Li S.
        • Garrett-Bakelman F.E.
        • Figueroa M.E.
        • Melnick A.
        • et al.
        methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles.
        Genome Biol. 2012; 13: R87
        • Amith M.T.
        • Fujimoto K.
        • Tao C.
        NET-EXPO: a gephi plugin towards social network analysis of network exposure for unipartite and bipartite graphs.
        HCI Int. 2019 Posters (2019). 2019; 1034: 3-12
        • Fornes O.
        • Castro-Mondragon J.A.
        • Khan A.
        • van der Lee R.
        • Zhang X.
        • Richmond P.A.
        • et al.
        Jaspar 2020: update of the open-access database of transcription factor binding profiles.
        Nucl. Acids Res. 2020; 48: D87-D92
        • Wisniewski J.R.
        • Zougman A.
        • Nagaraj N.
        • Mann M.
        Universal sample preparation method for proteome analysis.
        Nat. Met. 2009; 6: 359-362