In this study, RNA-Seq and bioinformatics analysis were performed on the mammary gland tissues of dairy goats in the LL, DP, and LG stages, and the gene expression profiles related to mammary gland development and lactation were mapped to construct mammary gland growth, proliferation, apoptosis, immunity, substance transport, synthesis, and metabolic-related networks. Prior to this study, studies have reported the transcriptome analysis of mammary gland tissues during the lactating and DP stages in pigs [18], sheep [25], and cattle [26]. The differences between this study and the above studies are that we selected three consecutive periods of mammary gland development. From the perspective of lactation, this is a process from the gradual termination of milk secretion to the absence of milk secretion and then to the restart of lactation. From the perspective of tissue structure, this is a process from mammary gland tissue involution to tissue remodeling, cell proliferation, and differentiation to form new mammary gland duct branches and acinar structure. Data analysis showed that a total of 1,381 DEGs were identified in the three mammary gland developmental stages, indicating that gene expression was indeed different among the three stages (Fig. 2). Perhaps the differential expression of these genes leads to differences in the morphology, structure and function of the mammary gland during the three stages. We identified the 20 DEGs with the highest expression in the mammary gland during the three stages, and these DEGs included CSN1S1 and CSN1S2, which encode alpha-casein, CSN2, which encodes beta-casein, CSN3, which encodes kappa-casein, and LALBA, which encodes alpha-lactalbumin. The expression levels of these gens were lower in the DP stage than in the LL and LG stages. The expression patterns of these genes found in this study were consistent with those of previous studies. [26, 27], and their high expression levels may be closely related to milk synthesis during the LL and LG stages. FGF2 and FGF11, as members of the fibroblast growth factor family, exhibited similar expression patterns, and their expression levels were lower in the LL stage than in the DP and LG stages. It has been reported that FGF2 promotes the proliferation of mammary epithelial cells, which is conducive to mammary gland remodeling and maintaining milk secretion [28]. RPS3A, RPS2, RPL19, RPLP0, RPL13A, RPL23, and RPS8 encode ribosomal proteins. Similarly, they also exhibited similar expression patterns, namely, expression levels during the LL stage were lower than those in the DP and LG stages. Therefore, we hypothesize that genes in the same family or genes with similar functions often have similar expression pattern in the mammary gland during different developmental stages. Other studies reported that in some cases, some gene family members had the same (or almost the same) sequence. These families allow large amounts of genes to be expressed in a short time as needed [29]. However, we used a limited number of goat mammary gland samples, and studies need to be conducted to verify this hypothesis.
The physiological structure of the goat mammary gland changes significantly during the three developmental stages. These changes are closely related to the specific expression of genes during different developmental stages [30]. According to the development/growth network (Fig. 5A), IGFBP5, a member of the insulin-like growth factor family, was associated with mammary gland development and morphogenesis and was highly expressed in the LL stage. It has been reported that IGFBP5 can accelerate the apoptosis of mammary epithelial cells and that the overexpression of IGFBP5 in vivo can lead to impaired mammary gland development [31]. Different from IGFBP5, IGFBP4 and IGFBP3 were expressed at low levels in the LL stage and at higher levels in the DP and LG stages. Our findings are consistent with those by Alla et al. in a study of the mammary gland tissues of adolescent, gestational, lactating and degenerative mice [32]. This again confirmed that insulin-like growth factor binding protein family members may have a potential regulatory relationship with mammary gland development. In addition, WNT5A and RTN4 genes involved in mammary epithelial cell proliferation had similar expression patterns, and all were expressed at low levels in the DP stage and high levels in the LG stage. It has been reported that WNT5A has an inhibitory effect on mammary gland development [33] and that RTN4 can inhibit hepatocellular carcinoma cell (HCC) proliferation and promote HCC apoptosis [34], which may counteract cell proliferation in LG stage. However, mitogen-activated protein kinase 1 (MAPK1), ESR1, cyclin D1 (CCND1), and mediator complex subunit 1 (MED1) were expressed at low levels in the LL stage and at high levels in the DP and LG stages. All these genes could promote the proliferation of mammary gland cells [35–38]. Therefore, RTN4 and WNT5A expression in the LG stage may be related to the regulation of mammary gland cell proliferation by other genes and to the differences in roles of these two genes in different tissues and cell types [35]. Transforming growth factor beta 3 (TGFB3) and transforming growth factor beta receptor 2 (TGFBR2) are both associated with the TGF-β signaling pathway, have become important molecular markers of breast cancer, and play important roles in the TGF-β pathway for mammary gland morphogenesis and lactation [39]. In addition, the elevated levels of prolactin and progesterone during pregnancy are conducive to mammary gland tissue remodeling and the formation of mammary gland lobe and duct branch structures [40]. However, we found that PGR and PRLR, as receptors for progesterone and prolactin, respectively, were highly expressed in mammary gland tissues during the LG stage, further demonstrating their promoting effects on mammary gland development, differentiation and lactation during the LG stage.
Apoptosis is essential to maintain the stability of the internal environment of the mammary gland, cell renewal and normal cell functions and activities [41]. According to the apoptosis network (Fig. 5B), members of the caspase family, CASP8, CASP9, and CASP3, are associated with extrinsic apoptotic signaling pathways, and their expression levels during the LL stage were much higher than those during the DP and LG stages. This is consistent with previous reports; that is, when the mammary gland degenerates, caspases are activated, and the expression levels of caspases are upregulated [42]. CASP8 and FADD like apoptosis regulator (CFLAR), galectin 3 (LGALS3), interferon alpha-inducible protein 6 (IFI6), lamin A/C (LMNA), janus kinase 2 (JAK2), TNF receptor superfamily member 1B (TNFRSF1B), secreted frizzled-related protein 1 (SFRP1), and protein tyrosine phosphatase receptor type C (PTPRC) are also associated with extrinsic apoptotic signaling pathways, and they are highly expressed in the DP stage. Among them, IFI6 encodes anti-apoptotic proteins that reside on the outer mitochondrial membrane, and interrupting the mitochondrial function of IFI6 may prevent breast cancer metastasis [43]. The constitutive activation of JAK2 in mammary epithelial cells can enhance signal transducer and activator of transcription 5 (Stat5) signal transduction, promote the generation of mammary acinus, inhibit apoptosis, and promote the develop of breast tumors [44]. SFRP1-encoded protein can block the Wnt signaling pathway, and low expression of SFRP1 can promote mammary epithelial cell proliferation and the branching of mammary ducts [45]. Based on existing reports and high expression levels in the DP stage, these three genes might be involved in the regulation of cell apoptosis in mammary gland tissue during remodeling in the DP stage. In addition, CFLAR-encoded proteins are regulators of cell apoptosis and are structurally similar to caspase-8. LGALS3 is a member of the galectin family, which encodes carbohydrate binding proteins, and members of this family have affinity for β- galactoside [46]. LMNA, TNFRSF1B and PPTRC are all associated with cell differentiation [47–49]. However, their role in the mammary gland is still unclear; therefore, further verification of their specific function in mammary gland development is needed. RACK1 is not only associated with the intrinsic apoptotic signaling pathway, but its expression level is in the top 20 during the DP and LG stages. The expression level of RACK1 in the LL stage was the lowest, while it was the highest in the LG stage. The scaffold proteins encoded by RACK1 can participate in the recruitment and assembly of a variety of signaling molecules, including protein kinase C [50], phosphodiesterases [51], Src tyrosine kinase [52], and insulin and insulin-like growth factors [53]. Therefore, RACK1 may participate in the apoptosis of mammary gland cells by regulating a variety of signaling molecules.
Immunomodulation is essential for mammary gland health and is the basis and guarantee for its normal function. On the one hand, the milk secreted by the mammary gland provides nutrition and immune protection for newborn pups; on the other hand, the mammary gland involution and DP stages are not only periods for mammary gland repair and renewal but are also critical periods for defending against various mammary diseases [54]. According to the immunity network (Fig. 5C), we found that genes associated with the regulation of complement activation, such as C5AR1, C7, C1S, SERPING1, CFB, C2, C1R, C1QB, C1QA, C1QC, and CFI, had similar expression patterns; they were all highly expressed in the DP stage. These results suggested that these genes might play a potential regulatory role in mammary gland immunity during the DP stage. C2, C1R, C1S, C1QA, C1QB, and C1QC are all associated with the classical immune pathway [55], have antibacterial functions, regulate immune adhesion, mediate inflammation, and play important roles in inherent immunity against microorganisms [56]. In addition, neutrophil-mediated immunity is also associated with innate immunity. lipocalin 2 (LCN2) was highly expressed in the LL stage. It is known that protein encoded by LCN2 is neutrophilic gelatinase-associated lipocalin, which inhibits the growth of bacteria through the chelation of iron-containing carriers, thereby playing a role in innate immunity [57]. In addition, C5AR1, S100A9, bone marrow stromal cell antigen 1 (BST1), fibrinogen like 2 (FGL2), PTPRC, cathepsin H (CTSH), and tubulin beta class I (TUBB) were highly expressed in the DP stage. C5AR1 encodes a receptor for the activated fragment C5a generated after the activation of complement and can be used as a chemokine to induce neutrophil phagocytosis to pathogens [58]. S100A9 encodes calcium and zinc binding proteins and plays an important role in the regulation of inflammatory processes and the immune response. It can induce neutrophil chemotaxis and adhesion and can enhance the anti-bactericidal activity of neutrophils through promoting phagocytosis by the activation of spleen associated tyrosine kinase (SYK), PI3K/AKT and extracellular signal-regulated kinase 1/2 (ERK1/2) [59]. BST1, FGL2, PTPRC, CTSH, and TUBB can also induce neutrophils to participate in the fight against pathogens, thus playing a potential protective role in mammary gland health during the DP stage. All three genes, peptidoglycan recognition protein 1 (PGLYRP1), beta-2-microglobulin (B2M), and LTF, were highly expressed in the LG stage. LTF (lactoferrin) is a transferrin protein that transfers iron to cells and controls the levels of free iron in blood and external secretions. It is widely present in the milk of humans and other mammals and has antibacterial activity. It is an important part of the nonspecific immune system [60]. In addition, humoral immune responses, B cell-mediated immunity, and T cell-mediated immunity also play important roles in the mammary gland. It is known that the protein encoded by JCHAIN can connect two monomeric units of IgM or IgA and induces a larger IgA polymer. This protein also helps to bind these immunoglobulins with secretory components to resist pathogens [61]. Both CD74 and CD81 are associated with B cell-mediated immunity, and CD74 is highly expressed in the DP stage, while CD81 is highly expressed in the LG stage. In summary, mammary gland immunity is associated with both the innate immune system and the specific immune system. The collaboration of these genes in the immune system during the three stages allows mammary gland immunity to proceed smoothly, thereby protecting mammary gland health.
The synthesis of proteins, lipids and lactose in milk begins with the transport of basic substances, such as amino acids, fatty acids and glucose, in plasma [62]. It can be seen from the transport network (Fig. 6A) that the members of the solute carrier family that were differentially expressed during the three stages were involved in fatty acid (SLC27A1, SLC27A3, and SLC27A4) amino acid (SLC1A5, SLC6A8, SLC6A14, and SLC25A29), and glucose (SLC2A4, SLC2A9, SLC2A10, and SLC2A8) transport. The members of these solute carrier families were highly expressed in the LL and LG stages, while low expression was observed in the DP stage. Compared with the gestation stage, the mRNA expression levels of SLC27A3 and SLC27A4 were upregulated 37-fold and 1.4-fold, respectively, during the lactation period [63]. In addition, SLC1A5, SLC2A4, and SLC2A8 were expressed at higher levels during the lactation stage than during the DP stage [26]. This indicates that the members of the solute carrier family are very important for substance transport to the mammary gland during the lactation stage, and their high expression during lactation is closely related to the higher lactation yield requirement. In addition, FABP3, PLIN2, MFSD2A, CD36 and other genes are also involved in fatty acid transport in the mammary gland and are highly expressed during the lactation period. By using proteomics, Zheng et al. found that the expression level of PLIN2 in the LL stage was lower than that in the peak period, which had a potential impact on milk fat synthesis in the mammary gland [64]. CD36 and FABP3 were significantly upregulated in the lactation period [65], and FABP3 is one of the most abundant isoforms in bovine mammary gland tissues; FABP3 expression is affected by the lactation period [66]. Proteins encoded by MFSD2A can mediate lipid transport, participate in the establishment of the blood-brain barrier and are required for brain growth and function [67]. In summary, these four genes might cooperate to participate in fatty acid uptake and transport in the mammary gland. In addition, we also found that ADP ribosylation factor like GTPase 6 interacting protein 1 (ARL6IP1) may be involved in amino acid transport; however, its expression pattern is different from that of solute carrier family members. ARL6IP1 can regulate amino acid transport and negatively regulate cell apoptosis by regulating the activity of caspase-9 (CASP9) to protect cell survival [68]. Therefore, ARL6IP1 is highly expressed in the LG stage and has low expression levels in the LL and DP stages. It may be involved in the transportation of amino acids in the LG stage; on the other hand, it may promote the proliferation and differentiation of mammary gland cells in the LG stage. The above data related to substance transport in the mammary gland indicate that during lactation, at different periods, the mammary gland has different demands for glucose, amino acids and fatty acids, and that the RNA expression of the related three types of transport proteins will also undergo proper adjustment and change, adapting to the lactation requirement of the mammary gland. The entire transport process is complex, which may require the coordination of the proteins encoded by multiple genes to achieve organized substance transport in the mammary gland.
The biosynthesis and metabolism of substances are the basic physiological processes of mammary gland development and lactation [69]. In addition to the genes encoding the main proteins in milk, such as casein (CSN1S1, CSN1S2, CSN2, and CSN3) and lactalbumin (LALBA), according to the biosynthesis network (Fig. 6B), a total of 12 genes (GLUL, OAT, CTH, BCAT2, CDO1, PSAT1, GPT2, PYCR1, glutaminase (GLS), phosphoglycerate dehydrogenase (PHGDH), phosphoserine phosphatase (PSPH), and asparagine synthetase (ASNS)) associated with cellular amino acid biosynthetic processes were identified. CDO1, PYCR1, PSAT1, BCAT2, ASNS and PSPH all had similar expression patterns and were highly expressed in the LL stage. CDO1 is known to be a key regulator of cysteine concentrations in cells and is identified as a candidate gene affecting milk protein traits [70]. PYCR1 encodes an enzyme that catalyzes the conversion of NAD(P)H-dependent pyrroline 5-carboxylate to proline [71]. In addition, proteins encoded by PSAT1 and PSPH are involved in serine synthesis [72]. BCAT2 encodes a branched-chain aminotransferase, and the encoded protein forms a dimer. This dimer catalyzes the first step in the production of branched-chain amino acids, such as leucine, isoleucine and valine. The abnormal expression of BCAT2 also affects cell proliferation and differentiation [73]. ASNS-encoded proteins are involved in the synthesis of asparagine [74]. The expression patterns of PYCR1, PSAT1, and ASNS were consistent with previous reports that the expression levels of these three genes during the lactation period are higher than those during the DP stage [26]. OAT, GPT2, PHGDH, and CTH are also involved in amino acid biosynthesis, and their expression levels were higher in the LG stage. These genes play important regulatory roles in amino acid biosynthesis in the mammary gland during the lactation period and provide enough amino acids in the mammary gland to achieve protein synthesis. The mammary gland is not only a place for protein synthesis but also an important place for lipid synthesis. It is known that the mammary gland of a lactating mouse synthesizes and secretes milk fat equivalent to its entire body weight in a 20-day lactation cycle, making the mammary gland one of the most active known lipid synthesis organs [62]. According to the biosynthesis network, we identified 12 genes associated with fatty acid synthesis. These include genes associated with fatty acid uptake in blood (LPL), the synthesis of long-chain fatty acids (TECR and LPGAT1) and short-chain fatty acids (ACSS1 and ACPS2), de novo fatty acid biosynthesis (ACACA and FASN), fatty acid desaturation (FADS1 and SCD), transcriptional regulation (INSIG1 and CD74), ketone body utilization (PDK4), and redox (GPX4). It has been reported that LPL, ACSS1, ACSS2, FASN, FADS1, and INSIG1 have similar expression patterns, all of which are highly expressed in the peak lactation stage, and the expression pattern of LPL in the bovine mammary gland is very similar to the lactation curve pattern, which may indicate that this gene plays an important role in maintaining milk synthesis [75]. SCD is a core enzyme in lipid metabolism and is the rate-limiting enzyme in the catalyzed synthesis of monounsaturated fatty acids from saturated fatty acid [76]. These lipid metabolism-related genes show different expression patterns during the three developmental stages, which are mutually adaptive to the fat synthesis function in different developmental stages. In addition, genes involved in carbohydrate synthesis and metabolism were also significantly enriched. GCK, ATF4, PGM1, and KAT2B are related to gluconeogenesis and glucose metabolic processes. GCK encodes a hexokinase family member, and hexokinase phosphorylates glucose to produce 6-phosphate glucose, which is the first step in most glucose metabolic pathways (Fig. 6C). The enzyme encoded by PGM1 is involved in the decomposition and synthesis of glucose. GCK and PGM1 are known to have higher expression levels during lactation than during the DP stage [77, 78]. The expression profiles of ATF4 and KAT2B in the mammary gland during different lactation stages were mapped for the first time. In summary, genes related to the biosynthesis and metabolism of basic substances such as glucose, amino acids and fatty acids in the mammary gland were screened using biosynthesis and metabolism networks, and their expression profiles in the mammary gland during three development stages were revealed. This study provides a reference for studying their roles in mammary gland development and lactation function in the future.
In the mammary gland, TFs interact with specific gene sequences (open or closed), and coordinated gene expression occurs at the correct time and in appropriate numbers. [79]. During mammary gland development and the lactation cycle, many site-specific TFs are essential for the proper growth, branching and involution of the mammary gland [80]. It is known that the overexpression of ATF4 inhibits the proliferation and differentiation of mammary epithelial cells, leading to lactation dysfunction and accelerated mammary gland involution [81]. Our results showed that ATF4 was highly expressed during the LL stage, possibly to promote the involution of mammary gland tissue and apoptosis during the LL stage (Fig. 7A-D). In addition, the expression of JUNB increases after forced weaning [82]. We also detected high JUNB expression in the DP stage, and we speculate that JUNB has effects on milk synthesis and secretion. CEBPG, FOSL2, and NFIL3 belong to the bZIP family, as do ATF4 and JUNB. bZIP TFs are present in all eukaryotes and are one of the largest families of dimeric TFs. It has been reported that this family of TFs regulates cell growth [83]. We found that this family of TFs had the greatest number of TFs with differential expression in the mammary gland; therefore, these TFs might be involved in mammary gland development and lactation regulation. In addition, four members of the bHLH family (ARNT, BHLHA15, BHLHE40, and MITF) were differentially expressed in the three developmental stages. ARNT is known to be necessary for mammary gland development in adolescent females [84]. BHLHE40 is associated with breast cancer [85]. However, the function of BHLHA15 and MITF in the mammary gland has not yet been reported; therefore, further study on their specific role in mammary gland development is warranted. In summary, these 35 TFs are differentially expressed in mammary gland tissues during the three developmental stages and may play potential regulatory roles in gene expression in the mammary gland, thereby affecting mammary gland development and the physiological process of lactation.