Transcriptome profiling of the nonlactating mammary glands of dairy goats reveals the molecular genetic mechanism of mammary cell remodeling

The mammary gland redevelops to the prepregnancy state during involution, which shows the mammary cells have the characteristics of remodeling. The rapidity and degree of mammary gland involution vary across species (e.g., between model organism mice and dairy livestock). However, the molecular genetic mechanism of involution and remodeling of goat mammary gland has not yet been clarified. This work investigated the structural changes and transcriptome characteristics of the mammary gland tissue of nonlactating dairy goats during the late lactation (LL), the dry period (DP), and late gestation (LG). Terminal deoxynucleotidyl transferase dUTP nick-end labeling (TUNEL) staining revealed significant changes in the structure of the nonlactating goat mammary gland, and obvious cell apoptosis occurred at LL and DP. Sequencing identified 1,381 genes that are differentially expressed in mammary gland tissue at the 3 developmental stages. Genes related to cell growth, apoptosis, immunity, nutrient transport, synthesis, and metabolism exhibited adaptive transcriptional changes to meet the needs of a new set of mammary gland lactation functions. The significant enrichment of Gene Ontology terms such as humoral immune response, complement activation, and neutrophil-mediated immunity indicates that the innate immune system plays an important role in maintaining the health of degenerative mammary glands and eliminating apoptotic cells. The peroxisome proliferator-activated receptor signaling pathway plays an important regulatory role in lipid metabolism, especially the adaptive changes in expression of genes encoded lipid transport and enzymes, which promote the formation of milk fat during the lactation. The mammary gland development gene module revealed that pregnancy hormone receptors, cell growth factors and their receptors, and genes encoding insulin-like growth factor binding proteins regulate the physiological process of mammary gland involution through adaptive transcriptional changes. Interestingly, ERBB4 was identified as the hub gene of the network that regulates mammary gland growth and development. Overexpression of ERBB4 in mammary epithelial cells cultured in vitro can reduce cell cycle arrest in G 1 /S phase and apoptosis by regulating the PI3K/Akt signaling pathway and promote the proliferation of mammary epithelial cells. The gene ERBB4 also affects the expression of genes that initiate mammary gland involution and promote mammary gland remodeling. These findings contribute to an in-depth understanding of the molecular mechanisms involved in mammary gland involution and remodeling.


INTRODUCTION
The mammary gland is an important exocrine gland in the breasts of mammals that is responsible for producing milk (Oftedal, 2002).Affected by hormonal changes during adolescence and pregnancy, the mammary gland undergoes many rounds of proliferation, expansion, and involution; thus, its development is unique (Brisken et al., 1999;Macias and Hinck, 2012;Inman et al., 2015).Mammary gland involution is a highly complex multistep process in which the lactating gland returns to a morphologically near-prepregnant state (Stein et al., 2007).There are 3 types of involution: gradual involution, senile involution, and active or initiated involution (Hurley, 1989).Mammary gland involution involves extensive tissue remodeling and is characterized by decreased synthesis of milk components, programmed cell death of massive secretory epithelial cells, increased permeability of tight junctions between alveolar epithelial cells, recruitment of immune cells, removal of apoptotic cells and stagnant milk, and refilling the mammary gland with differentiated fat cells (Nonnecke and Smith, 1984;Hurley, 1989;Strange et al., 1992).

Transcriptome profiling of the nonlactating mammary glands of dairy goats reveals the molecular genetic mechanism of mammary cell remodeling
The molecular mechanisms of mammary gland involution, involving signaling molecules and pathways, have been widely discussed (Hurley, 1989;Watson and Kreuzaler, 2011;Jena et al., 2019;Zhao et al., 2019).Rodents are the most widely studied animals, and their mammary gland involution is divided into 2 stages (Watson, 2006).These 2 stages are controlled by the gradual increase in death signals and the loss of survival factors.The lysosomal-mediated cell death (LCD) pathway and the mitochondrial-mediated intrinsic apoptotic pathway are the main apoptotic pathways at the first and second stages of mammary gland involution, respectively (Lund et al., 1996;Clarkson et al., 2004).The first stage of involution is controlled by signals from the local mammary gland (Li et al., 1997).After cessation of milking, stagnant milk in the mammary gland triggers an increase in the expression of death-inducing signaling molecules through complex mechanisms (Stein et al., 2007).For example, the BCL2 Associated X, apoptosis regulator (BAX) (Rucker et al., 2011), tumor necrosis factor (TNF) (Kojima et al., 1996), signal transducer and activator of transcription 3 (STAT3) (Hughes and Watson, 2018), transforming growth factor β-3 (TGFB3) (Flanders and Wakefield, 2009), Interleukin 10 (IL10) (Hwang et al., 2016), p53 (Jerry et al., 1999), leukemia inhibitory factor (LIF) (Schere-Levy et al., 2003) were significantly upregulated at the beginning of involution.Milk stasis also induces the loss of STAT5A and STAT5B phosphorylation, thereby destroying the main pathway of prolactin signaling (Tan and Nevalainen, 2008).The second stage is characterized by tissue remodeling caused by proteases and hormonal changes at the systemic level (Li et al., 1997).The systems and pathways involved in the second stage of involution mainly include plasminplasminogen-plasminogen activator, metalloprotease, and prolactin (Lund et al., 1996).In addition, the mammary gland involution also involves the immune system (Atabai et al., 2007), IGF-IGF binding protein system (Neuenschwander et al., 1996), PI3K/Akt signaling pathway (Schwertfeger et al., 2001;Rädler et al., 2017), the connection between the mammary epithelial cells and the matrix (Nguyen and Neville, 1998), Zinc transporter (SLC30A2; Hennigar et al., 2015), and other molecular mechanisms.
Transcriptome sequencing provides a high-throughput method for studying the molecular mechanisms of mammary gland development.Currently, a large number of studies have been carried out on transcriptome analysis of mammary gland or milk, involving goats (Crisà et al., 2016;Shi et al., 2016;Li et al., 2020), sheep (Paten et al., 2015;Wang et al., 2020), cattle (Bionaz et al., 2012;Zheng et al., 2018;Bhat et al., 2019;Fan et al., 2021), and pig (Shu et al., 2012;Palombo et al., 2018) at different stages of mammary gland development.Researches on the expression profile of mammary gland tissue after cessation of milking provide important clues to the molecules that may mediate epithelial cell death and tissue remodeling during involution (Clarkson et al., 2004;Stein et al., 2004;Shi et al., 2016).However, the degree and rapidity of mammary gland involution differ between different species, especially between ruminants and monogastric animals (Zhao et al., 2019).It has been reported that the involution is caused by the increase in pressure in the goat's mammary glands (Fleet and Peaker, 1978), whereas in rodents, the stop of nursing stimulation seems more critical.Unlike rodents, the lobular alveolar structure in mammary gland is maintained during the involution of dairy cows.Upon cessation of milking, adipocytes extensively regenerate and fill the mammary gland of mouse (Elias et al., 1973;Hovey and Aimo, 2010).The involution of adipose tissue in the mammary gland is different in dairy cows from mice.Dairy cow mammary epithelium experienced limited degeneration during the involution.The udders of cows at this time were rumored to resemble "wrung-out dishrag," reflecting the mammary glands were not to recover significant lipid deposits during involution (Hovey and Aimo, 2010).However, compared with the mice and the dairy cows, studies related to the transcriptome characteristics of the nonlactating mammary gland and the molecular genetic mechanism of mammary gland involution in goats are still insufficient.Dairy goats provide the milk for humans.They have important economic value and are ideal models for studying mammary gland development.In an earlier transcriptomics analysis of the mammary glands of goats during the early, peak, and late stages of lactation, we identified factors that regulate mammary gland development and the lactation cycle and determined their expression patterns (Ji et al., 2019(Ji et al., , 2020)).However, these findings represent only the tip of the iceberg in understanding the complex physiological processes that occur during mammary gland development and lactation and emphasize the need for further exploration.Therefore, in the current study, transcriptome analysis of the mammary gland tissues of nonlactating goats was conducted as a means of investigating the molecular genetic mechanism of mammary gland involution and remodeling.

Ethics Statement
The methods and experimental procedures used in this study were approved and directed by the animal protection and ethics committee of Shandong Agricul-  .figshare .14973861 .v1;Xuan, 2021a).All experimenters were required to undergo rigorous training and self-protection to minimize the pain experienced by the experimental animals.

Experimental Animal Treatment and Sample Collection
This research used healthy, disease-free Laoshan dairy goats raised at the Qingdao Aote goat breeding farm (Qingdao, Shandong, China).The goats were provided with the farm's standard feeding and management environment.All goats were of third parity, and the range of age was 3 yr, 8 mo to 4 yr, 1 mo; the range of goat BW was 50.16 to 56.41 kg (specific information on the goats is available in Table S1 in the previous study by Xuan et al., 2020a).Mammary gland tissues were collected during late lactation (LL; n = 3, 240 d postpartum), the dry period (DP; n = 3, 300 d postpartum), and late gestation (LG; n = 3, 140 d after mating).The specific tissue collection methods were as follows.After intravenous injection of pentobarbital sodium (100 mg/ kg), the muscles relaxed, and the heart and respiration were arrested; the goats were then quickly dissected, and mammary gland tissues were harvested.The harvested tissues were immediately placed in RNase-free cryotubes and stored in liquid nitrogen.

Differential Expression Analysis
First, genes with low expression were filtered based on the fragments per kilobase of transcript per million mapped reads (FPKM).Before differential expression analysis, principal component analysis and intersample correlation analysis of gene expression levels were performed.Differential expression analysis was then performed using DESeq2 software (Love et al., 2014).Genes were considered differentially expressed between 2 groups if the false discovery rate (FDR) was <0.05 and the log2-fold change was ≥1.The differentially expressed genes (DEG) between different comparison groups (LL vs. DP, DP vs. LG, LL vs. LG) were analyzed, and differential gene expression was displayed using volcano plots and heat maps.A Venn plot was used to show the distribution of DEG between different comparison groups in each developmental period and in different developmental periods.The DEG were then ranked according to their expression levels at each developmental stage.The top 20 DEG for each developmental period (LL, DP, and LG) were identified and are displayed using histograms.

Gene Ontology Annotation and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
The DEG identified in each group were annotated using the Gene Ontology (GO) database (Ashburner et al., 2000) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000), and specific enrichment analysis was realized by the R package clusterProfiler (Yu et al., 2012).In addition, to interpret more clearly whether the enriched GO terms or signaling pathways are overall upregulated or inhibited during mammary involution, we run the same analyses with up-and downregulated DEG in each comparative group.Goat genes annotated in the GO and KEGG databases were used as background genes for enrichment analysis.The GO terms and KEGG pathways that met FDR <0.05 were selected for downstream analysis.To further clarify and refine the functions of genes, we first classified GO terms and KEGG pathways into multiple categories, including apoptosis-related, biosynthesis-related, growth and development-related, metabolism-related, substance transport-related, and other important function-related terms.The GO terms and KEGG pathways related to mammary gland development and lactation function are presented as bar graphs and bubble plots, respectively.

Functional Classification and Analysis of Expression Patterns of DEG
To further clarify the relationship between specific genes and important GO terms, we constructed geneconcept networks, including development/growth, apoptosis, immunity, transport, biosynthesis, and metabolism networks, using the R package clusterProfiler.The patterns of gene expression changes in the network are shown in heat maps.The FPKM values of the gene expression levels were normalized based on the Z-score before the heat maps were plotted, and the genes were then hierarchically clustered.The expression heat map was drawn by the R package pheatmap, and the R package TCseq was used to draw a line graph of gene expression patterns.

Protein-Protein Interaction Network Analysis and Identification of Hub Genes Among DEG Related to Mammary Gland Development
To further screen for hub genes related to mammary gland development, protein interactions for which the combined score was ≥0.4 were screened for the genes in the development/growth network using STRING software (version 11.5;Szklarczyk et al., 2019).The network was drawn using Cytoscape software (Shannon et al., 2003).Each note's degree was calculated using Cytoscape's built-in software NetworkAnalyzer, and genes of the appropriate degree were defined as hub genes.In addition, KEGG pathway enrichment analysis was performed on genes in the protein-protein interaction (PPI) network using the R package cluster-Profiler, and the pathways and genes with FDR <0.05 were screened for display.

Screening and Analysis of Transcription Factors
Transcription factors (TF) are key regulators of gene expression (Latchman, 1997).To classify the TF found among the DEG, we used the AnimalTFDB 3.0 website to annotate the DEG (Hu et al., 2019).The transcription factors information of Capra hircus provided on the website AnimalTFDB 3.0 (http: / / bioinfo .life.hust.edu.cn/static/ AnimalTFDB3/ download/ Capra _hircus _TF) were used as the annotation background.The numbers of TF belonging to different families were then counted, and the expression patterns of these differentially expressed TF were analyzed.

Hematoxylin and Eosin Staining and Terminal dUTP Nick-End Labeling Assay
Mammary gland tissue obtained at each of the 3 developmental stages (LL, DP, LG) was rinsed with 1× PBS (pH = 7.4), and tissue blocks (1 × 1 × 0.3 cm) were prepared using surgical scissors.Four percent paraformaldehyde solution (Beyotime) was used to fix and preserve the tissue blocks.The fixed mammary gland tissue was rinsed in running water for 24 h, and the fixation fluid was removed.Pieces of tissue of appropriate size were cut and dehydrated by sequential incubation in 50, 70, 80, 90, and 100% alcohol.Xylene solution was used to render the dehydrated mammary tissue transparent, and finally, the tissue block was embedded in paraffin.The paraffin-embedded mammary gland tissue was sectioned, and routine hematoxylin and eosin staining of the sections was performed.In addition, 20 μg/mL DNase-free proteinase K (Beyotime) was added dropwise to other mammary gland tissue sections, and those sections were incubated at 37°C for 15 min.After washing out the proteinase K with 1× PBS (pH = 7.4), 50 μL of terminal dUTP nick-end labeling (TUNEL) detection solution (Beyotime) was added to each sample and the sections were incubated at 37°C in the dark for 40 min.The cells were washed 3 times with 1× PBS (pH = 7.4) for 5 min each time.The slices were then incubated in 4′,6-diamidino-2-phenylindole (DAPI) staining solution (Beyotime) for 10 min.Images were captured using Pannoramic 250 Flash Series digital scanners (3DHISTECH); Ca-seViewer (3DHISTECH) and ImageJ software (version 1.53;Schneider et al., 2012) were used to analyze the images.Six microscopic fields (20× magnification) in each stained section were randomly selected, and the cells in each field were counted.The ratio of TUNELlabeled cells to DAPI-counted nuclei was defined as the mammary epithelial cell apoptosis rate.

Goat Mammary Epithelial Cell Culture
Goat mammary gland tissue was surgically obtained.After the mammary gland tissue was rinsed with sterile PBS, adipose and connective tissue were removed using ophthalmic scissors, and the tissue was cut into fragments (0.5-1 mm 3 ).Goat primary mammary epithelial cells were then cultured at 37°C in 5% CO 2 using the tissue block method (Wang et al., 2013).Gibco Dulbecco's modified Eagle's medium was used for cell culture.Nutrient mixture F-12 (DMEM/F-12) medium was supplemented with 10% fetal bovine serum (Gibco, Thermo Fisher Scientific), 0.25 mM hydrocortisone, 5 mg/mL insulin (Sigma), 50 U/mL penicillin/streptomycin, and 10 ng/mL epidermal growth factor 1 (Invi-Xuan et al.: TRANSCRIPTOME PROFILING OF NONLACTATING GOAT MAMMARY GLANDS trogen).Subsequently, mammary gland epithelial cells were seeded in 6-well culture plates at a seeding density of 3×10 5 .When the cells reached approximately 75% confluence, they were prepared for transfection.

Cell Proliferation Assay
Goat mammary epithelial cells were seeded in 96well plates at a density of 8,000 cells per well.When the cells had reached 75% confluence, cell transfection with the overexpression vector and the interference vector was performed.The CCK8 method was used to measure cell proliferation.The specific steps included harvesting of the cells at 0, 24, 48, and 72 h, addition of 10 μL of CCK8 solution (Beyotime) to each well, and incubation of the cells at 37°C for 1 h.The absorbance values of the treated wells were recorded at 450 nm using a SpectraMax ID3 microplate reader (Molecular Devices).GraphPad Prism Version 8 software was used to create a line chart showing the cell proliferation activity.

Apoptosis Assay
Mammary epithelial cells were seeded in 6-well culture plates at a density of 3×10 5 cells/well.Forty-eight hours after transfection, 5 μL of annexin V-fluorescien isothiocyanate (FITC) and 10 μL of propidium iodide staining solution (Invitrogen) were added to each well, and the plates were incubated for 10 min at room temperature in the dark.A BD Flow Cytometer (BD Biosciences) was then used for apoptosis detection.FlowJo 7.6.1 software (BD Biosciences) was used to analyze the flow cytometry results.Western blotting was used to detect the protein expression levels of caspase-3, Bax, and BCL-2 in cells 48 h after transfection.

Western Blotting and Immunohistochemistry of Mammary Gland Tissue
Radioimmunoprecipitation assay (RIPA) lysis solution (Beyotime) was used to extract total protein from mammary gland tissue and from epithelial cells.A bicinchoninic acid (BCA) kit (Beyotime) was used to measure the protein concentration of the extract at 562 nm on a microplate reader (Molecular Devices).A 10% separating gel and a 5% concentrating gel were prepared (Beyotime).Twenty micrograms of total protein from each sample was separated by SDS-PAGE, and the separated proteins were transferred to polyvinylidene fluoride (PVDF) membranes.Each PVDF membrane was blocked in buffer containing 5% skim milk powder for 1.5 h and incubated at 4°C for 12 h in buffer containing a primary antibody against the target protein.After 3 washes with 1× Tris-buffered saline containing Tween 20, the PVDF membrane was incubated with the appropriate secondary antibody for 50 min at room temperature.The membrane was then infiltrated with ECL chemiluminescence buffer (Beyotime) for 1 min, and the immunoreactive protein was detected using an Azure 300 Chemiluminescent Western Blot Imaging System (Azure Biosystems).ImageJ 1.48 software (National Institutes of Health) was used to calculate the protein expression level, and GAPDH was used as an internal reference for expression level correction.Immunohistochemistry staining was also performed on mammary gland tissue slice specimens obtained at the 3 developmental stages.The specific steps were as follows: the mammary gland slices were baked in a thermostat at 50°C for 20 min and then immersed in xylene solution for 25 min to remove paraffin.The slices were then hydrated through immersion in a graded alcohol series.After washing out the alcohol with PBS, the slices were incubated with 3% H 2 O 2 for 10 min to inactivate endogenous peroxidase activity and boiled for 10 min in 0.01 M sodium citrate buffer (pH = 6.0).After natural cooling, blocking solution was used for 20 min.The slices were then incubated with a primary antibody against ErbB4 at 4°C for 12 h.A negative control was generated by replacing the primary antibody with 1× PBS solution.The sections were washed 3 times in 1× PBS at room temperature for 7 min each time and incubated with a goat antirabbit IgG secondary antibody at 37°C for 1 h.Streptavidinperoxidase was then added dropwise, and the sections were incubated at room temperature for 30 min, rinsed 3 times with PBS for 15 min each time, and stained with 3,3′-diaminobenzidine for 5 min.Pannoramic 250 Flash Series digital scanners (3DHISTECH) were used to capture images of the sections, and CaseViewer software (version 2.4; 3DHISTECH) was used to analyze the images (20× field of view).The antibodies used in this study included anti-ErbB4 (1:1,500 dilution, ab219208, Abcam), anti-CCND1 (1:1,500 dilution, #55506, CST), anti-GAPDH (1:2,000 dilution, #2118, CST), anti-PI3K (1:1,000 dilution, #4257, CST), anti-p-PI3K (1:1,000 dilution, #17366, CST), anti-Akt (1:1,000 dilution, #4691, CST), anti-p-Akt (1:1,000 dilution, #4060, CST), anti-caspase-3 (1:2,000 dilution, #14220, CST), anti-BCL-2 (1:1,000 dilution, ab182858, Abcam), anti-Bax (1:1,000 dilution, #2772, CST) and HRP-conjugated goat antirabbit IgG (1:2,000, CW0103, CWBIO).

Real-Time Quantitative Reverse-Transcription PCR Assay
To verify the accuracy and reliability of the transcriptome data, we randomly selected 15 DEGs for real-time quantitative reverse-transcription PCR (RTqPCR) verification.In addition, to determine the effect of ErbB4 overexpression on mammary gland involution, we selected 11 genes related to mammary gland involution for quantitative measurement of RNA transcripts.The National Center for Biotechnology Information Primer-BLAST was used to design gene primer sequences (Ye et al., 2012).The RNA was extracted using the TRIzol method.Reverse transcription and PCR amplification were performed using a One Step TB Green PrimeScript RT-PCR Kit (Takara).The RT-qPCR was performed on a LightCycler 96 realtime thermal cycler (Roche).Use NormFinder software to select the most stable gene pair from the 6 genes as the reference genes for RTqPCR assay.Calculate the primer amplification efficiency by establishing a standard curve.The geometric mean of 2 carefully selected reference genes (GAPDH and MRPL39) was used as an accurate normalization factor (Vandesompele et al., 2002).Relative gene expression calculated by the improved Pfaffl Method (Pfaffl, 2001; specific calculation steps at https: / / toptipbio .com/qpcr -multiple -reference -genes/ ).The primer sequences, efficiency, and stability value of reference genes are all provided in Supplemental Table S1 (https: / / doi .org/ 10 .6084/m9 .figshare.14967357.v5;Xuan, 2021q).

Statistical Analysis
Student's t-test was used to test the significance of differences between the experimental group and the control group; P < 0.05 was taken to indicate a sig-nificant difference.Tukey's honestly significant difference test was used to analyze gene expression levels at different developmental stages.Unless otherwise noted, the plots presented in this work were generated using R software (Version 3.6.0;https: / / www .r-project .org/).

Basic Statistical Analysis of Sequencing Results
Each library obtained an average of 153 million clean reads through transcriptome sequencing, and the alignment rate with the reference genome was ≥95% (Table 1).All data have been uploaded to the Gene Expression Omnibus database (https: / / ncbi .nlm.nih.gov/geo/ ) under the accession number GSE185981.The principal component analysis and cluster analysis results between samples showed that 3 samples clustered together in the same period (Supplemental Figure S2, https: / / doi .org/ 10 .6084/m9 .figshare.14966640.v1;Xuan, 2021b).These results indicate that the sequencing data are stable and reliable, and the differences between the repeat samples from the same period are small.There were differences between samples from different periods.Through quantitative analysis, 14,523 nonredundant expressed genes were identified in the mammary gland (FPKM ≥0.5).The numbers of genes expressed at the DP and LG stages were 50 and 44% higher, respectively, than the number of genes expressed at the LL stage (Table 1).The results of differential expression analysis showed a total of 1,381 genes whose expression levels changed significantly during the 3 developmental stages (Supplemental Table S2, https: / / doi .org/ 10 .6084/m9 .figshare.14967381.v1;Xuan, 2021r).Among them, the number of DEG found in the comparison between DP and LG was the largest (Figure 1A).In addition, 62 genes were differentially expressed in the comparison among the 3 developmental stages (Figure 1B).The global expression heat map (Figure 1C) and the volcano map (Supplemental Figure S3, https: / / doi .org/ 10 .6084/m9 .figshare.14966715.v1;Xuan, 2021c) clearly show that the expression of 1,381 genes changed apparently during the 3 periods.Based on gene expression level at each mammary gland developmental stage (LL, DP, LG), the top 20 DEG were screened.The expression of genes encoding milk components (CSN1, CSN2, CSN1S1, CSN1S2, LALBA, BTN1A1), growth factors (FGF2, FGF11), enzymes (MGAT2), serum amyloid A protein (LOC106503728, LOC100860781, LOC102168428) and ribosomal proteins (RPS2, RPS3A, RPS8, RPL13A, RPL19, RPL23, RPLP0) changed apparently (Figure 2D-F).In addition, the use of RTqPCR to verify 15 DEG (Supplemental Figure S4, https: / / doi .org/ 10 .6084/m9 .figshare.14966787.v2,Xuan, 2021d) showed that the differential expression log2-fold change in expression was highly correlated with the RNA-seq data (correlation value = 0.83, P-value <0.001), indicating that the RTqPCR data are consistent with the RNA-seq data.These results indicate that major changes occur at the gene transcriptional level during the involution and remodeling of mammary gland tissue from the LL to the LG stage, reflecting the stage specificity of gene expression in mammary gland tissue.

Transcription Factor Analysis
A total of 35 TF were identified among a total of 1,381 DEG (Supplemental Figure S5, https: / / doi .org/ 10 .6084/m9 .figshare.14966817.v1;Xuan, 2021e).The TF were divided into 16 families, but family information for 4 of the TF was not clear.The basic leucine zipper (bZIP) family had the largest number of TF; it included the 5 TF ATF4, CCAAT enhancer binding protein gamma (CEBPG), FOS like 2 (FOSL2), JunB proto-oncogene, AP-1 transcription factor subunit (JUNB), nuclear factor, and interleukin 3 regulated (NFIL3).The expression heat maps indicated that the expression patterns of CEBPG and ATF4 were similar; the expression levels of those TF were higher at the LL stage than at the DP and LG stages.The expression patterns of NFIL3, JUNB and FOSL2 were similar, and their expression levels were higher at the DP stage than at the LL and LG stages.Four TF in the basic helixloop-helix (bHLH) family, such as aryl hydrocarbon receptor nuclear translocator (ARNT), A15 (BHLHA15), BHLHE40, and microphthalmia-associated transcription factor (MITF), were found in the TF analysis.The genes MITF, BHLHE40, and ARNT showed similar expression patterns, and the expression levels of all 4 identified bHLH family TF were highest at the DP stage.Estrogen receptor 1 (ESR1), progesterone receptor (PGR), and estrogen-related receptor α (ESRRA) all belong to the ESR-like family.These receptors also had similar expression patterns, with high expression in the LG stage and relatively low expression at the LL stage.Similar expression patterns of most TF in the same family were also found for TF in the homeobox and zf-C2H2 families.Please refer to Supplemental Table S3 (https: / / doi .org/ 10 .6084/m9 .figshare.14967402.v1;Xuan, 2021s) for specific information on the TF analysis.These results indicate that the expression of TF changes in the 3 stages; the observed changes may play a regulatory role in gene expression during the process of mammary gland involution and remodeling.

GO and KEGG Pathway Enrichment Analysis
We provide a global overview of the functions of 1,381 DEG at different developmental stages and the possible signaling pathways involved.The DEG were classified into 6 functional categories, including apoptosis, growth and development, immune function, substance transport, metabolism, and biosynthesis, as well as other GO terms (Supplemental Figure S6, https: / / doi .org/ 10 .6084/m9 .figshare.14966841.v1;Xuan, 2021f), and with respect to KEGG signaling pathways (Supplemental Figure S7A, https: / / doi .org/ 10 .6084/m9 .figshare.14966859.v2;Xuan, 2021g).According to the functional enrichment analysis of DEG in the LL versus DP group (Figure S7B: downregulated DEG in LL vs. DP), upregulated genes during the DP are enriched to the biological process of cell adhesion and aging, such as negative regulation of cell adhesion, cell substrate adhesion, regulation of cell substrate adhesion, extracellular matrix organization, and aging.The KEGG pathway analysis shows that upregulated genes during the DP are involved in multiple immune-and disease-related signaling pathways, such as staphylococcus aureus infection, complement and coagulation cascades, and antigen processing and presentation.Gene function analysis of DEG in DP vs. LG shows (Figure S7C: upregulated DEG in DP vs. LG): Upregulated genes in the DP involve tissue remodeling and external encapsulating structure organization.Upregulated genes in the LG (Figure S7C: downregulated DEG in DP vs. LG) are involved in material metabolism (fatty acid metabolic process, galactose metabolism), material transport (long-chain fatty acid import across plasma membrane, organic anion transport), mammary gland LG (Supplemental Figure S7D), it was found that upregulated genes at LG are closely related to the proliferation of mammary epithelial cells (epithelial cell proliferation, response to epidermal growth factor, and PI3K-Akt signaling pathway).Many signaling pathways related to tumorigenesis and cancer are also enriched, suggesting that mammary cells proliferate actively at the LG.In addition, upregulated genes at the LL are enriched in apoptosis-related signal-ing pathways.Overall, at the transcriptional level, the DEG found for the nonlactation period may be involved in the regulation of a variety of physiological activities, showing the diversity involved in the physiological process of mammary gland development.At the same time, this finding implies that a variety of physiological processes must be coordinated during mammary gland development.Please refer to the specific GO enrichment results (Supplemental  LG, LL vs. LG).These genes are then ranked according to the expression level at each stage, and the top 20 genes at each stage are selected to display.FPKM = fragments per kilobase of transcript per million mapped reads.The data are presented as the mean ± SEM, n = 3.

Structural Characteristics of Mammary Gland Tissue and Cell Apoptosis
The hematoxylin and eosin staining results showed (Figure 2A) that at the LL stage, the mammary acinus was atrophied with a reduced volume and irregular shape, and the acinus cavity contained flocculent protein secretion products.At the DP stage, the mammary acinus is replaced by extensive connective tissue.The epithelial cells of the mammary gland at the LG are round and full, surrounded by a regular acinus that is filled with milk and lipid droplets, and the acinus structure is clear and complete.The TUNEL results showed that extensive apoptosis occurred in mammary tissue cells at the LL and DP stages, and there was a similar degree of apoptosis and a higher rate of apoptosis in these 2 periods than at the LG stage (41.2 ± 6.5-fold and 31.8 ± 3.1-fold, respectively; Figure 2B).The sequencing results (Supplemental Figure S8, https: / / doi .org/ 10 .6084/m9 .figshare.14966880.v2;Xuan, 2021g) showed a total of 39 apoptosis-related genes.Twentyeight genes were significantly enriched in the intrinsic apoptotic signaling pathway, and 11 genes were significantly enriched in the extrinsic apoptotic signaling pathway (FDR < 0.05).In addition, 7 proapoptotic genes with high expression at the LL stage, 6 proapoptotic genes with high expression at the DP stage, and 8 antiapoptotic genes with high expression in the LG stage were found (Figure 2C).The PPI network identified 7 hub genes.The hub genes CASP3 and ATF4 were highly expressed in LL, whereas P4HB was highly expressed in LG.These results indicate that extensive apoptosis occurs at the LL and DP stages.Apoptosis, which helps remove aging and damaged mammary cells, is an important step in mammary gland involution and remodeling.

Analysis of Immune Function Genes Expressed in the Mammary Gland
The GO enrichment results (the top 20 most significantly enriched terms) showed that the enrichment significance and the number of genes enriched were higher for the category "humoral immune response" than for the category "lymphocyte-mediated immunity" (Supplemental Figure S9, https: / / doi .org/ 10 .6084/m9 .figshare.14967174.v2;Xuan, 2021h).Moreover, among the top 20 enrichment terms, GO terms related to innate immunity were 3 times more highly enriched than those related to adaptive immunity.All 10 hub genes identified by the PPI network were related to innate immunity.This finding indicates that both innate and adaptive immunity are involved in protecting the health of the mammary gland during the process of involution and remodeling.Nevertheless, the innate immune response plays a more obvious role and is more active than the adaptive immune response.In addition, this study found that a total of 14 genes were enriched in the category "complement activation," and 13 of those genes were highly expressed at the DP stage.This result suggests that the complement system plays an important role in protecting mammary gland health at the DP stage.In addition, 12 genes were significantly enriched in GO terms associated with neutrophils.Four and 8 of these genes were highly expressed at the LL and LG stages, respectively.Higher rates of mammary cell apoptosis were detected during these 2 periods, and neutrophils are known to have strong chemotactic and phagocytic functions.Therefore, we speculate that these genes may play an important role in stimulating neutrophils to clear apoptotic mammary cells.In summary, we showed the characteristics of the expression patterns of 33 immune-related genes at 3 developmental stages of the mammary gland.Genes related to the innate immune system play an important role in mammary gland involution and remodeling; the action of their gene products is conducive to the formation of a first line of defense for the mammary gland against pathogen invasion.These genes also play an important role in activating the adaptive immune system and eliminating apoptotic cells.

Analysis of Genes Related to Substance Transport, Synthesis, and Metabolism
Based on the GO annotation results, this study classified 209 DEG with respect to 3 cellular components: 97 genes related to protein metabolism, 50 genes related to carbohydrate metabolism, and 79 genes related to lipid metabolism (Supplemental Figure S10, https: / / doi .org/ 10 .6084/m9 .figshare.16977448.v3;Xuan, 2021i).In the analysis of lipid metabolism-related genes, the peroxisome proliferator-activated receptor (PPAR) signaling pathway was significantly enriched, and the number of enriched genes was the largest (Supplemental Figure S11, https: / / doi .org/ 10 .6084/m9 .figshare.14974989.v3;Xuan, 2021j).The PPI network shows that the hub genes SCD, ACSL1, and LPL are involved in the PPAR signaling pathway, and they are all highly expressed at the LL and LG stages.Through the expression profile of the 12 genes in the PPAR signaling pathway, it can be seen that the expression of the genes related to fatty acid transport (SLC27A1, SLC27A4, SLC27A5, FABP3) and the enzymes that catalyze fat metabolism reactions (ACSL1, HMGCS1, LPL, PCK2, SCD, SCD5) at the DP are lower than that at the LL and LG stages (Supplemental Figure S12, https: / / doi .org/ 10 .6084/m9 .figshare.16975639.v1;Xuan, 2021k).However, the PPAR family members and their co-activators have not detected significant changes at their expression levels in the 3 periods.Carbohydrate-related pathways include glycolysis/gluconeogenesis, carbon metabolism, galactose metabolism, the pentose phosphate pathway, and other pathways (Supplemental Figure S13, https: / / doi .org/ 10 .6084/m9 .figshare.16985740.v3;Xuan, 2021l).Among them, 3 hub genes (GCK, PGM1, UGP2) are related to the regulation of galactose metabolism.Amino acids and proteins are the building blocks of life,

Xuan et al.: TRANSCRIPTOME PROFILING OF NONLACTATING GOAT MAMMARY GLANDS
the physiological processes related to AA and protein metabolism were significantly enriched (Supplemental Figure S14, https: / / doi .org/ 10 .6084/m9 .figshare.14967252;Xuan, 2021m).This study also demonstrated the expression patterns of genes during mammary gland involution in GO terms related to post-translational protein modification, amino acid transmembrane transport, protein folding.These results indicate that the mammary gland faces great challenges related to nutrient acquisition and utilization during the transition from dry to lactation.The flexible changes in the expression of metabolism-related genes at the transcriptional level meet the needs for metabolism homeostasis of material and energy in mammary gland tissue.
According to the expression level of gene, the metabolic genes are divided into 3 clusters (Supplemental Figure S15, https: / / doi .org/ 10 .6084/m9 .figshare.14973984.v3;Xuan, 2021n).The genes in cluster 1 are expressed at a higher level at the LL, whereas the genes in cluster 2 and cluster 3 are expressed at higher levels during the DP and LG, respectively.In addition, the number of genes in cluster 1 and cluster 3 is 1.8-fold and 2.2-fold of that in cluster 2, respectively (Supplemental Table S6, https: / / doi .org/ 10 .6084/m9 .figshare.16974619.v2;Xuan, 2021v).This difference reflects the decrease in substance metabolism at the DP stage, a decrease that may be related to the involution of mammary gland tissue, the slow metabolism of substances and energy, and the cessation of lactation.Material transport also account for an important part of the acquisition and processing of major nutrients.This finding shows that members of the solute carrier protein family contribute to the transport of nutrients in the mammary gland and that the genes in this family also display adaptive changes in expression in response to changes in physiological functions (Supplemental Figure S16, https: / / doi .org/ 10 .6084/m9 .figshare.14967297.v2;Xuan, 2021o).

Analysis of Genes Related to Mammary Gland Development and Cell Growth
A total of 53 DEG were related to mammary gland development and growth (Figure 3).Among them, hormones play an important role in the regulation of mammary gland development and involution.The prolactin signaling pathway, growth hormone synthesis, secretion and action, thyroid hormone signaling pathway, and other physiological processes associated with hormones were significantly enriched (Figure 3B-E).The genes PRLR, ESRRA, ESR1, and PGR were highly expressed in the LG (Supplemental Figure S17A, https: / / doi .org/ 10 .6084/m9 .figshare.14967324.v2;Xuan, 2021p); this enrichment may be related to the establishment and maintenance of pregnancy.This study found that MMP2 and MMP14 were highly expressed at the DP stage (Supplemental Figure S17B).These results indicate that the matrix metalloproteinase system is activated during the DP stage and that it induces proteolysis of extracellular matrix components, leading to extensive apoptosis and degeneration of mammary cells.The expression levels of cell growth factors, cell growth factor binding proteins, and cell growth factor receptors also changed.The PPI network identified 7 hub genes (Figure 3A), of which PIK3CA, AKT1, FN1, FGF2, CCND1, ERBB4, and MAPK1 are all related to the PI3K/Akt signaling pathway (Figure 3C).The ERBB signaling plays an important role in tissue development and homeostasis.Increasing evidence has shown that ERBB4 regulates the elongation of mammary ducts, but this process requires cell proliferation and epithelial cell expansion (Jones et al., 1999;Muraoka-Cook et al., 2008).We found that at the transcription level, the ERBB4 gene and other PI3K/Akt signaling pathwayrelated genes are highly expressed at LG (Figure 3F).However, the function and mechanism of ERBB4 in controlling mammary epithelial cell function are still unclear.Therefore, the role of the ERBB4 gene in mammary cells in vitro was explored in depth.

Expression of ErbB4 in Goat Mammary Gland
The expression of ErbB4 at the protein level was consistent with its expression at the transcriptional level, and it was highly expressed at both levels at the LG stage (Figure 4).This expression pattern was contrary to the number of apoptotic cells and the expression of proapoptotic genes at the different stages of mammary development.This finding indicates that ErbB4 expression may be related to mammary cell proliferation and apoptosis.The immunohistochemistry results showed that ErbB4 was localized in the nucleus and cytoplasm of mammary epithelial cells at the LL and LG stages.At the DP stage, ErbB4 was not clearly localized in the nucleus and was mainly expressed in the cytoplasm.This finding indicates that the location of the protein encoded by ERBB4 differs in mammary gland tissue cells at different developmental stages.

Overexpression of ErbB4 in Vitro Promotes Proliferation of Mammary Epithelial Cells by Activating the PI3K/Akt Signaling Pathway
To verify the effect of ErbB4 on the proliferation of mammary epithelial cells, we interfered with or overexpressed ErbB4 protein in cultured mammary epithelial cells in vitro and measured cell viability at different time points by CCK8 assay.The western blotting results showed that ErbB4 was expressed at a signifi-cantly lower level (55.6 ± 15.6% lower; P < 0.05) in the interference group than in the control group and at a significantly higher level (67.8 ± 18.6% higher; P < 0.05) in the overexpression group than in the control group (Figure 5A).The CCK8 results showed obvious differences in cell proliferation activity at 48 and 72 h after cell transfection (Figure 5C).At 48 h, compared with the control group, cell proliferation activity in the interference group was significantly reduced by 30.7 ± 4.5%, whereas in the overexpression group, it was significantly increased by 65.8 ± 15.8% (P < 0.05).In addition, the western blotting results showed that interference with ErbB4 decreased the p-PI3K/PI3K and p-Akt/Akt ratios compared with those observed in the control group, whereas overexpression of ErbB4 increased the p-PI3K/PI3K and p-Akt/Akt ratios.These data suggest that in vitro overexpression of ErbB4 promotes the proliferation of mammary epithelial cells through activation of the PI3K/Akt signaling pathway.

Overexpression of ERBB4 Inhibits Apoptosis of Mammary Epithelial Cells and Reduces Cell Cycle Arrest at the G 1 /S Phase
The flow cytometry results (Figure 6A) showed that at 48 h after transfection, the apoptosis rate in the ErbB4 interference group was 38.6% ± 8.4% higher than that in the control group.Conversely, the apoptosis rate in the overexpression group was 56.9 ± 7.3% lower than that in the control group.The western blotting results (Figure 6B) showed that the Bcl-2/Bax ratio was significantly lower in the interference group than in the control group (P < 0.01).In contrast, the Bcl-2/Bax ratio was significantly higher in the overexpression group than in the control group (P < 0.05).Overexpression of ErbB4 can also effectively reduce the expression of caspase 3. Cell cycle analysis (Figure 6C) showed that, compared with the control group, interference with ErbB4 increased the number of cells in G 1 /S phase, decreased the number of cells in the G 2 /M phase and reduced the expression of CCND1.In contrast, overexpression of ErbB4 increased the expression of CCND1 and caused more cells to enter the G 2 /M phase.In general, overexpression of ErbB4 activates the PI3K/Akt signaling pathway, inhibits the apoptosis of mammary epithelial cells, changes the cell cycle by upregulating CCND1, and reduces the number of cells arrested in G 1 /S phase and is thus conducive to the proliferation of mammary epithelial cells.

Effect of ErbB4 on the Genes that Initiate Mammary Gland Involution
Compared with the control group (Figure 7), interference with ErbB4 significantly increased the transcriptome levels of SATAT3, BCL2L11, CAPN2, and LGALS3 in mammary epithelial cells cultured in vitro, whereas it significantly decreased the expression of STAT1 and STAT5B.Stat signaling is known to play an important role in the initiation of mammary gland involution (Hennighausen et al., 1997).Therefore, it is speculated that ErbB4 may be involved in the process of inducing mammary gland involution.In addition, interference with or overexpression of ErbB4 in mammary epithelial cells had no effect on MFGE8 or ITGA6, members of the matrix metalloproteinase family (MMP2, MMP14), or their inhibitors (TIMP1).

DISCUSSION
This study mainly reported the transcriptome characteristics of nonlactating mammary gland and the molecular genetic mechanism of mammary gland involution in dairy goats.By transcriptome sequencing of goat mammary gland tissue, 1,381 DEG were identified at 3 mammary gland developmental stages (Figure 1).(A) Analysis of apoptotic goat mammary epithelial cells.Flow cytometry was used to detect the apoptosis of mammary epithelial cells after overexpression and interference with ErbB4.Forty-eight hours after the cells were transfected, the cells were stained with annexin V-FITC/ propidium iodide (PI), incubated in the dark for 10 min, and then detected using a BD flow cytometer.The bar graph shows the statistics of apoptotic cells in the different groups (si-ErbB4 vs. si-NC, ErbB4 vs. NC).(B) Analysis of apoptosis-related protein expression.The expression of proapoptotic (caspase-3, Bax), antiapoptotic (Bcl-2), and cyclin (CCND1) genes was detected by Western blotting.The bar graphs show the ratio of Bcl-2/Bax and the expression of Caspase-3 and CCND1, respectively.(C) Analysis of cell cycle assays.Forty-eight hours after cell transfection, the cells were stained with PI and incubated in the dark for 10 min for cell cycle detection using BD flow cytometry.The bar graph shows the proportion of cells at different stages of the cell cycle.si-NC: native control of small interfering RNA (empty vector of siRNA); si-ErbB4: the small interfering RNA targeting ErbB4; NC: native control (empty vector of pcDNA3.1);ErbB4: ErbB4 overexpression vector.*P < 0.05; **P < 0.01; data are presented as the mean ± SEM (n = 3).This finding indicates that there are significant changes in gene expression at the transcriptional level during mammary gland involution.Among these genes, some of those with high expression are related to protein synthesis (Figure 1D), and others were found to be related to mammary gland growth.The expression of genes encoding milk proteins (CSN1S1, CSN1S2, CSN2, CSN3, LALBA) clearly decreased during the DP stage because lactation did not occur during this time.This is consistent with previous studies on the transcriptome analysis of the mammary gland tissue of dairy cows during the lactation and DP (Dado-Senn et al., 2018;Dai et al., 2018;Lin et al., 2019).Similarly, genes encoding ribosomal subunits and responsible for AA transport are low expressed in the mammary gland tissue of goats during the DP (Figure 1, Supplemental Figure S14).In addition, 68 DEG that were upregulated at the LG stage were also closely related to the ribosome (Supplemental Figure S7C: downregulated DEG in DP vs. LG).As we all know, the ribosome is responsible for completing the translation process from RNA to protein in the central dogma of molecular biology.Our research found that the ribosome-related genes are not actively expressed during the DP, but are highly expressed at the LG stage (Supplemental Figure S7E).This indicates that, on the one hand, it may be related to the synthesis of more milk protein during lactation; on the other hand, this expression characteristic may be related to rapid cell proliferation and differentiation occurring in the mammary gland at LG, which requires the synthesis of more tissue structural proteins.Interestingly, studies of rabbit and sheep mammary glands support the first hypothesis.Because measurements of ribosomal proteins in rabbit and sheep mammary tissues show that the number of polyribosomes is increased in lactation mammary tissues compared with nonlactation mammary tissues (Denamur, 1974).However, the decrease of protein synthesis machinery in the mammary gland during lactation is a common feature between mice and cows.This is explained by the mammary gland preferentially encoding the mRNA of proteins related to milk synthesis and secretion at the translation level, rather than the mRNA of non-milk-specific proteins (Bionaz et al., 2012a).Therefore, it is an adaptation of the gene that encodes the ribosome to make better milk proteins.At the same time, it also indicates that the action efficiency of ribosomes in the mammary gland is altered at different lactation stages of the same species and in various species at the same lactation stage.
Protein folding and modification are essential for converting newly synthesized proteins into biologically functional forms (Braakman and Hebert, 2013).This study showed the transcriptional level expression characteristics of genes related to protein folding and modification during the involution of mammary glands (Supplemental Figure S14).The high expression of coding chaperone genes (e.g., DNAJC3, HSPA5) during lactation is related to the promotion of protein translation (Paten et al., 2015).Prolyl 4-hydroxylase subunit beta (P4HB), also known as protein disulfide bond isomerase, is a multifunctional protein that catalyzes the formation and rearrangement of disulfide bonds (Xiong et al., 2018).It can be used as a molecular chaperone to improve misfolded proteins in response to endoplasmic reticulum stress.Endoplasmic reticulum oxidoreductase 1 alpha (ERO1A) is a hypoxia-induced endoplasmic reticulum oxidase that can be regulated by the translation and folding of oxidized proteins.In addition, the high expression of ERO1A at the LL and DP may also be involved in the regulation of cell apoptosis (Li et al., 2009a).The protein encoded by PRDX4 is an antioxidant enzyme.The high expression of PRDX4 during lactation plays a role in regulating mammary cell protection and preventing oxidative stress.In addition, this study found genes related to protein modification, which are essential for enriching protein functions.At the same time, genes related to the unfolded protein response (UPR) were also enriched.The accumulation of milk protein in secretory mammary epithelium may generate stress signals that trigger involution, which may be through the integration of UPR signals, autophagy, and apoptosis during mammary gland involution (Wärri et al., 2018).The above illustrates the coordination of many aspects from the transportation of amino acids to protein translation and folding and protein modification to meet the needs of protein synthesis changes in the process of mammary gland involution.
In the mammary gland, TF interact with specific gene sequences (open or closed), and coordinated gene expression occurs at the correct times and at the appropriate levels (Mitchell and Tjian, 1989).The STAT are actively involved in various stages of mammary gland development; 7 STAT family members (STAT1, STAT2, STAT3, STAT4, STAT5A, STAT5B, and STAT6) have been identified to date in mammalian cells (Darnell et al., 1994).STAT3 is activated at the first stage of mammary gland involution in mice, and it promotes the LCD pathway (Li et al., 1997).We found that STAT1 was highly expressed only during the DP and LG stages (Supplemental Figure S5).This result differs from that of a previous report indicating that the expression of pSTAT1 is highest in adolescent glands and during degeneration and low during pregnancy and lactation.In addition, many recent reports indicate that STAT1 can affect mammary tumors, and the mechanism of its action in mammary development remains unclear (Haricharan and Li, 2014).The differences in expression patterns of the STAT family are influenced by the species differences between goats and mice and indicate the diversity of molecular functions of STAT.In addition, we found that the transcription levels of bZIP family members and bHLH family members changed during mammary gland involution.We found that ATF4 was highly expressed at the LG stage; this high expression is likely closely related to the fact that ATF4 overexpression has been reported to inhibit the proliferation and differentiation of mammary cells and accelerate mammary gland involution (Yonekura et al., 2018).These results enrich our knowledge of the regulatory role of TF in mammary gland development.
An important characteristic of mammary gland involution is the recruitment of immune cells (Farhadian et al., 2018;Boutinaud et al., 2019;Zhao et al., 2019).Immune cells (such as neutrophils and macrophages) in cow's milk will rise rapidly in the early stage of involution and remain until delivery (Ollier et al., 2013).At the same time, the immune signals are also changing during the involution process (Bu et al., 2017;Wang et al., 2020), the sequencing results could support this (Supplemental Figure S9).The complement system plays an important role in host defense mechanisms such as immune lysozyme action, immune adhesion, immune agglutination, and enhanced phagocytosis (Loor et al., 2011;Merle et al., 2015).Concentrations of complement are lowest in milk secreted by healthy mammary glands during lactation and highest in mammary secretions obtained during LL and the DP and in colostrum (Rainard, 2003;Ezzat Alnakip et al., 2014, Katsafadou et al., 2019).In this study, high expression of genes related to the complement system was detected in degenerative mammary gland tissue, a finding that also supports the above statement.Because the mammary gland involution process requires the removal of large amounts of apoptotic cells and milk components, it is reasonable to hypothesize that specialized phagocytes of the innate immune system play a crucial role in this process (Paape et al., 2000).We found that neutrophilrelated genes were highly expressed at the DP stage.The neutrophil-specific chemokine GRO1/CXCL1 was upregulated 5-fold on the first day of mammary involution in mice.The upregulation of the neutrophil marker LRG was highest on the second day after forced weaning (Stein et al., 2004), suggesting that neutrophils and their related genes experience adaptive changes in the process of mammary gland involution.On the one hand, this finding can help clear aging and apoptotic mammary gland cells; on the other hand, it can help the animal resist disease and protect the health of the mammary gland.This is consistent with the previously reported features and functions of immune cells during mammary gland involution of goats (Loor et al., 2011;Lin et al., 2015), sheep (Tatarczuch et al., 2002;Wang et al., 2020), and bovine (Dado-Senn et al., 2018).
In the face of the challenges related to nutrient uptake and metabolism during mammary gland involution, genes related to the metabolism of mammary substances have undergone adaptive changes.On the one hand, such changes are reflected in the changes in the expression of genes related to the biosynthesis of basic nutrients and metabolism at the transcriptional level (Supplemental Figure S10); on the other hand, members of the solute carrier protein family make important contributions to the transport of amino acids, fatty acids, and carbohydrates (Supplemental Figure S16).In addition, the PPAR signaling pathway plays an important regulatory role in milk fat synthesis and metabolism of mammary epithelial cells (Bionaz et al., 2013;Shi et al., 2013;Liu et al., 2016).The expression characteristics and effects of the PPAR isotypes are different in ruminants and monogastric animals (Bionaz et al., 2013).Among them, PPARγ (gene symbol PPARG) plays an important role in milk fat synthesis (Shi et al., 2013).Bovine mammary gland tissues have higher PPARG expression during pregnancy and lactation (Bionaz and Loor, 2008).The expression of PPARG was also detected in goat mammary glands, albeit at a significantly lower level compared with the cattle (Bernard et al., 2013).The expression of PPARG in the mammary glands of mice decreases from pregnancy to lactation (Rudolph et al., 2007), whereas the expression of PPARG in the mammary glands of pigs is not affected by the lactation cycle (Shu et al., 2012).Our research found that the PPAR family members (PPARA, PPARD, PPARG) had no significant differences in transcription levels during the mammary gland involution.On the one hand, this may be related to the characteristics of adipose tissue changes during goat mammary gland involution.On the other hand, it indicates that the expression of PPAR family members may be regulated by other factors, such as prolactin receptor signaling (Viengchareun et al., 2008) and noncoding RNA (Portius et al., 2017).However, the genes encoding lipid transport proteins and lipases have undergone significant changes during the mammary gland involution.They are highly expressed at the LL and LG, and decreased at the DP.This indicates that the PPAR signaling pathway might affect the lipid synthesis and metabolism during involution by regulating the transport of circulating or stored lipids in the body and coordinating various enzyme-catalyzed reactions (Gbaguidi et al., 2002;Shi et al., 2013;Zhang et al., 2021).In addition, the expression levels of carbohydrate-and protein-related genes were low at the DP stage, which is possibly related to the decrease in energy metabolism that occurs during the DP stage.In conclusion, the changes in the transcriptome levels of metabolism-related genes at different developmental stages reflect the adaptability of mammary gland metabolism.Knowledge of these changes also provides a reference for the reasonable supply of animal nutrients in the process of mammary gland involution.
Apoptosis is an important step in the process of mammary gland involution and remodeling, but there are differences among species (Watson, 2006;Stein et al., 2007).In terms of tissue structure, the involution of the mammary gland of dairy cows does not result in the collapse and death of a large number of acinar cells.The acinus structure presents heterogeneity; most cells maintain contact with the basement membrane, and all are closely connected with neighboring cells (Capuco and Akers, 1999).This connection may be related to the overlap between the gestation period and the DP of dairy cows.It is affected by pregnancy hormones and local signals that decrease the rate of apoptosis of mammary gland cells.In the physiological process, the number of apoptotic cells and the degree of apoptosis are significantly lower in bovine mammary gland cells than in the mammary gland cells of mice (Zhao et al., 2019).However, our study found that the goat mammary gland experienced significant apoptosis at both the LL and DP stages (Figure 2).This finding was also supported by the high expression of apoptosis-related genes at the LL and DP.Two distinct processes have characterized mammary gland involution in mice (Watson, 2006).However, we did not continuously sequence the mammary gland tissue during involution, which places some limitations on the ability to provide a detailed description of the characteristics of gene transcriptome changes during mammary gland involution.
Genes and important pathways involved in mammary gland development are altered during involution.We identified 8 hub genes through a PPI network and found that ErbB4 plays an important role in regulating the physiological process of mammary gland involution through the PI3K/Akt signaling pathway (Figure 3).The receptor tyrosine-protein kinase ErbB4 is required for STAT5A activation in the mammary gland during pregnancy (Clark et al., 2005).We found high expression of ErbB4 at both the mRNA and protein levels at the LG stage (Figure 4).The ErbB receptor tyrosine kinase family is usually associated with increased epithelial cell growth and malignant transformation and progression (Hardy et al., 2010).In contrast, ErbB4/HER4 exhibits unique properties.Proteolysis is performed in 2 steps, releasing an 80 kDa nuclear-localized tyrosine kinase that acts as a signal transduction mechanism, slows down growth, and stimulates mammary cell differentiation (Muraoka-Cook et al., 2008).The PI3K/ Akt pathway is a key signaling node for the expansion and differentiation of luminal mammary epithelium (Wickenden and Watson, 2010).Sequencing results showed that ErbB4 was involved in the regulation of the PI3K/Akt signaling pathway.To further confirm the function of ErbB4 in mammary cell remodeling, we overexpressed the ErbB4 gene in goat mammary epithelial cells.We found that ErbB4 can promote the proliferation of mammary epithelial cells by activating the PI3K/Akt signaling pathway (Figure 5).In addition, activation of Akt is involved in the regulation of the cell cycle and proliferation and anti-apoptosis (LaRocca et al., 2011;Xuan et al., 2020b).On overexpression of ErbB4 in cultured mammary epithelial cells in vitro, the Bcl-2/Bax balance was disrupted, the number of apoptotic cells was reduced, and G 1 /S phase arrest of the cell cycle was slowed (Figure 6).These results indicate that ErbB4 is a favorable factor for promoting mammary gland remodeling.
Regulation of the process of mammary gland involution involves multiple factors.However, there are differences in mammary involution between ruminants and rodents.In this study, the influence of ErbB4 overexpression on the expression of representative genes in the first and second stages of goat mammary involution was studied by referring to previous mouse mammary involution-related genes (Jena et al., 2019).Interference with ErbB4 in mammary epithelial cells cultured in vitro can trigger the high expression of STAT3, the initiation gene of mammary involution (Figure 7).The involution of the mammary gland is characterized by extensive apoptosis of epithelial cells and activation of STAT3 (Haricharan and Li, 2014).The genes STAT3 and STAT5 have a mutual activation pattern throughout the mammary gland development cycle, which indicates that STAT5 may be a survival factor of MECs and that STAT3 may be a death factor of MECs.We found that ErbB4 overexpression can activate the expression of STAT5A and STAT5B.These results further suggest that ErbB4 is important for the initiation of mammary gland involution.Interestingly, the gain or loss of ErbB4 functionality does not affect matrix metalloproteinase members, suggesting that ErbB4 plays a significant role in regulating intramammary signaling pathways but has little effect on exocrine signaling (MMP3, MMP9, TIMMP1) (Muraoka-Cook et al., 2008).
This study has some limitations.One is the problem of statistical power.Three animals were used as replicates at each stage of mammary gland development.At least 5 transcriptome replicates are recommended for better testing power (Yu et al., 2017).Second, the differences in genetic information between individual animals may interfere with the results of differential expression analysis in the mammary gland at different periods (Conesa et al., 2016;Hekman et al., 2016).In addition, this study evaluated only 3 time points during mammary gland development (long time intervals).Although differences in the function and structure of mammary glands were found, mammary gland involution and remodeling are continuous processes, and the use of more sampling and sequencing time points might help explain the molecular change mechanism of this process and might better reflect the continuous characteristics of dynamic gene changes during mammary gland involution.

CONCLUSIONS
In this study, we analyzed the mammary gland tissue structure and transcriptome characteristics of nonlactating goats during the mammary gland involution.Cell apoptosis during the LL and the DP was essential for removing mammary senescent cells and maintaining the stability of the mammary gland internal environment.Gene expression patterns related to cell growth, apoptosis, immunity, material transport, biosynthesis, and metabolism were significantly altered during mammary gland involution.The innate immune response is more active than adaptive immunity, and it plays an important role in maintaining mammary gland health during the involution period.Compared with other nutrients, genes related to amino acid transport, ribosomal subunits, protein folding, and modification are more active at the LG and lactation.The PPAR signaling pathway is an important pathway for regulating lipid metabolism and production of milk fat.In addition, ErbB4 is conducive to mammary gland tissue remodeling by regulating the PI3K/Akt signaling pathway, which affects the expression of genes that initiate mammary gland involution.These findings provide new insights into the regulation of cell fate, polarity, branching morphogenesis, and functional organ involution during mammary gland development.
of China (31672401), National Key Research and Development Plan of China (2018YFD0501906), Shandong Provincial Natural Science Foundation of China (ZR2014CM029), the Shandong Provincial Modern Agriculture Industry Technology System (SDAIT-10-01), and funds of the Shandong "Double Tops" Program (SYL2017YSTD12).The authors have not stated any conflicts of interest.

Figure 1 .
Figure 1.Results of gene differential expression analysis.(A) Numbers of differentially expressed genes (DEG) in each group.LL = late lactation; DP = dry period; LG = late gestation.(B) Venn diagram of DEG.(C) Heatmap of all DEG.(D) Top 20 DEG at each developmental stage of the mammary gland.The genes are DEG obtained from pairwise comparisons (LL vs. DP, DP vs.LG, LL vs. LG).These genes are then ranked according to the expression level at each stage, and the top 20 genes at each stage are selected to display.FPKM = fragments per kilobase of transcript per million mapped reads.The data are presented as the mean ± SEM, n = 3.

Figure 3 .
Figure 3. Analysis of mammary gland growth and development functional module.(A) Protein-protein interaction (PPI) network analysis of differentially expressed genes (DEG) related to mammary gland development.(B) Dot plot of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis results.(C) Networks between DEG and KEGG pathways.(D) Networks between DEG and Gene Ontology (GO) terms.(E) Bar plot of GO annotation results.(F) Expression patterns of DEGs related to the PI3K/Akt signaling pathway at 3 different developmental stages of the mammary gland.LL = late lactation; DP = the dry period; LG = late gestation; FPKM = fragments per kilobase of transcript per million mapped reads.*P < 0.05; **P < 0.01; the data are presented as the mean ± SEM (n = 3).

Figure 4 .
Figure 4. Expression of ErbB4 and CCND1 in nonlactating mammary glands.(A) Western blot analysis of the expression of ERBB4 and CCND1; GAPDH was used as an internal reference for relative quantitative analysis.(B) Detection of ErbB4 expression in mammary gland tissues at 3 stages by immunohistochemistry. NC: negative control (1× PBS instead of primary antibody); ErbB4: mammary tissue sections were treated with ErbB4 primary antibody.LL = late lactation; DP = the dry period; LG = late gestation.*P < 0.05; **P < 0.01; the data are presented as the mean ± SEM (n = 3).

Figure 5 .
Figure 5.Effect of overexpression or interference with the expression of ErbB4 in mammary epithelial cells on the PI3K/Akt signaling pathway.(A) Western blotting was used to detect the protein expression of ERBB4 after overexpression or interference in mammary epithelial cells.(B) The effect of overexpression or interference of ERBB4 on the PI3K/Akt signaling pathway was detected by western blotting.(C) The proliferation of mammary epithelial cells was detected by CCK8 assay.si-NC: native control of small interfering RNA (empty vector of siRNA); si-ErbB4: small interfering RNA targeting ErbB4; NC: native control (empty vector of pcDNA3.1);ErbB4: ErbB4 overexpression vector.*P < 0.05; **P < 0.01; the data are presented as the mean ± SEM (n = 3).

Figure 6 .
Figure 6.Effect of overexpression or interference with the expression of ErbB4 in mammary epithelial cells on apoptosis and the cell cycle.(A)Analysis of apoptotic goat mammary epithelial cells.Flow cytometry was used to detect the apoptosis of mammary epithelial cells after overexpression and interference with ErbB4.Forty-eight hours after the cells were transfected, the cells were stained with annexin V-FITC/ propidium iodide (PI), incubated in the dark for 10 min, and then detected using a BD flow cytometer.The bar graph shows the statistics of apoptotic cells in the different groups (si-ErbB4 vs. si-NC, ErbB4 vs. NC).(B) Analysis of apoptosis-related protein expression.The expression of proapoptotic (caspase-3, Bax), antiapoptotic (Bcl-2), and cyclin (CCND1) genes was detected by Western blotting.The bar graphs show the ratio of Bcl-2/Bax and the expression of Caspase-3 and CCND1, respectively.(C) Analysis of cell cycle assays.Forty-eight hours after cell transfection, the cells were stained with PI and incubated in the dark for 10 min for cell cycle detection using BD flow cytometry.The bar graph shows the proportion of cells at different stages of the cell cycle.si-NC: native control of small interfering RNA (empty vector of siRNA); si-ErbB4: the small interfering RNA targeting ErbB4; NC: native control (empty vector of pcDNA3.1);ErbB4: ErbB4 overexpression vector.*P < 0.05; **P < 0.01; data are presented as the mean ± SEM (n = 3).

Figure 7 .
Figure 7. Effect of ERBB4 on genes related to mammary gland involution.Reverse transcription quantitative PCR was used to detect the expression of mammary gland involution-related genes in mammary epithelial cells after overexpression or interference with ERBB4.The bar graphs show the relative expression of 11 genes at transcription level in different groups (si-ErbB4 vs. si-NC; ErbB4 vs. NC).si-NC: native control of small interfering RNA (empty vector of siRNA); si-ErbB4: the small interfering RNA targeting ErbB4; NC: native control (empty vector of pcDNA3.1);ErbB4: ErbB4 overexpression vector.* indicates P < 0.05, ** indicates P < 0.01; data are presented as mean ± SEM (n = 3).

Table 1 .
Xuan et al.: TRANSCRIPTOME PROFILING OF NONLACTATING GOAT MAMMARY GLANDS Basic statistical analysis of the sequencing data development (mammary gland development, prolactin signaling pathway), ribosomal function (ribosome biogenesis, ribosome), and other physiological processes and signaling pathways.According to functional analysis of DEG in LL vs.