Abstract

Background

   Lactogenesis includes two stages. Stage I begins a few weeks before
   parturition. Stage II is initiated around the time of parturition and
   extends for several days afterwards.

Methodology/Principal Findings

   To better understand the molecular events underlying these changes,
   genome-wide gene expression profiling was conducted using digital gene
   expression (DGE) on bovine mammary tissue at three time points (on
   approximately day 35 before parturition (−35 d), day 7 before
   parturition (−7 d) and day 3 after parturition (+3 d)). Approximately
   6.2 million (M), 5.8 million (M) and 6.1 million (M) 21-nt cDNA tags
   were sequenced in the three cDNA libraries (−35 d, −7 d and +3 d),
   respectively. After aligning to the reference sequences, the three cDNA
   libraries included 8,662, 8,363 and 8,359 genes, respectively. With a
   fold change cutoff criteria of ≥2 or ≤−2 and a false discovery rate
   (FDR) of ≤0.001, a total of 812 genes were significantly differentially
   expressed at −7 d compared with −35 d (stage I). Gene ontology analysis
   showed that those significantly differentially expressed genes were
   mainly associated with cell cycle, lipid metabolism, immune response
   and biological adhesion. A total of 1,189 genes were significantly
   differentially expressed at +3 d compared with −7 d (stage II), and
   these genes were mainly associated with the immune response and cell
   cycle. Moreover, there were 1,672 genes significantly differentially
   expressed at +3 d compared with −35 d. Gene ontology analysis showed
   that the main differentially expressed genes were those associated with
   metabolic processes.

Conclusions

   The results suggest that the mammary gland begins to lactate not only
   by a gain of function but also by a broad suppression of function to
   effectively push most of the cell's resources towards lactation.

Introduction

   Lactating cows are generally dried off simply by stopping the milking
   process approximately two months before the next parturition. The
   mammary gland then undergoes an involution process, which is marked by
   the cessation of secretory activity and the reabsorption of milk
   residue, followed by a relatively static period. The dry period has
   been proven important for dairy cows. Omitting or shortening the dry
   period imposes negative effects on mammary health and milk yield in the
   next lactation [31][1]–[32][3]. The mammary gland does not resume its
   activity until approximately 2 to 3 weeks before the next parturition,
   when dramatic changes occur to prepare for profuse milk secretion after
   parturition. Lactogenesis is defined as the process from the resumption
   of mammary activity until profuse milk secretion and is divided into
   two stages [33][4]. Stage I is the period from the resumption of
   mammary activity to the time of parturition and is characterized by
   mammary differentiation, proliferation, and progressive expression of
   milk protein, as well as the secretion of pre-colostrum. Stage II is
   initiated around the time of parturition and extends for several days
   afterwards. This stage is characterized by the closure of the tight
   junctions between alveolar cells and the formation and secretion of
   colostrum and milk. The mammary gland is the only organ that
   experiences regular proliferation and involution cycles after maturity,
   which makes it an ideal model for the study of organ development.
   Knowledge of the molecular events driving lactogenesis in dairy cows
   has contributed not only to the understanding of organ development but
   also to the development of new technologies in the management and
   breeding of dairy cattle. To date, knowledge about this aspect of the
   mammary gland has mainly come from studies using mammary cell lines and
   genetically modified mice [34][5]–[35][7]. Some proteins involved in
   lactogenesis and their related signaling or metabolic pathways have
   been identified [36][8], [37][9]. It has been suggested that there is
   no sudden transcriptional switch around the time of parturition.
   Preparation of the gland for lactation includes modifications to the
   transcriptional program, but the onset of lactation appears to be
   primarily controlled by post-transcriptional mechanisms [38][10]. Lemay
   et al [39][11] analyzed the microarray data sets of mammary gland RNA
   samples collected from FVB mice at 10 time points during mammary
   development, and the results indicated SAM68 (an RNA-binding
   transduction protein and a putative regulator of mRNA splicing,
   transduction and nuclear export) to be an important
   post-transcriptional regulator of both milk secretion and mammary cell
   survival during lactation. Finucane et al [40][12] studied the
   molecular events in stage II of lactogenesis using Affymetrix GeneChip
   Bovine Genome Array, and found that genes associated with cell cycle
   and proliferation were downregulated. It was consistent with the result
   of Sorensen et al [41][13], which indicated that the mammary glands of
   cows proliferated mainly in late pregnancy and almost ceased
   proliferating after parturition. However, in dairy goats, mammary
   growth continued into early lactation, peaking at day 5 of lactation
   [42][14]. In addition, Suchyta et al [43][15] compared the microarray
   generated transcript profiles of liver, spleen, thymus, adrenal, ileum,
   and lymph tissues collected from a 3 month old Holstein steer and of
   mammary tissues collected from pre-pubertal and post-pubertal Holstein
   heifers, and identified a putative set of 16 genes being preferentially
   expressed in the developing mammary gland. In stage II of lactogenesis,
   the mammary gland undergoes a set of developmental processes that lead
   to the secretion of colostrum, and then milk. The major compositions of
   milk are lactose, protein and fat. Bionaz and Loor [44][9], [45][16]
   evaluated the expression of 44 genes involved in bovine mammary milk
   protein synthesis and 45 genes involved in milk fat synthesis via
   quantitative PCR. The results of those studies supported a pivotal role
   for the concerted action of PPARG, PPARGC1A, and INSIG1 in the
   regulation of milk fat synthesis and a central role of amino acid and
   glucose transporters and insulin signaling through mTOR in the
   regulation of protein synthesis in the bovine mammary gland.

   Currently, little is known concerning the molecular events underlying
   lactogenesis. Comparisons of the transcriptomes of pre- and
   post-parturition mammary glands of dairy cows by Kiera et al [46][12],
   provided only some insights into the molecular events of stage II. To
   gain information on what occurs during stage I, the transcriptomes of
   the relatively static mammary gland and that of stage I should be
   compared. Therefore, mammary gland samples from Holstein cows at
   approximately −35 d, −7 d and +3 d relative to parturition were
   collected for the present study. Deep sequencing of the transcriptomes
   in these samples was performed using the newly developed digital gene
   expression (DGE) method, which was thought to be more reliable,
   repeatable, and precise compared to previous microarray technology
   [47][17].

Materials and Methods

Ethics Statement

   Our study had been approved by the Institutional Animal Care and Use
   Committee of Shandong Agricultural University and performed in
   accordance with the “Guidelines for Experimental Animals” of the
   Ministry of Science and Technology (Beijing, China).

Mammary tissue collection

   Mammary gland samples were collected from 18 Holstein cows reared in
   dairy farm No. 1 of the Jiabao dairy company, Jinan, China. They were
   all in the second parity with the same age. There were 6 cows in each
   group. Mammary gland tissue samples were collected from cows at three
   different time-points wherein distinct differences exist in the
   internal status of the gland: one, 35±2 days before parturition (−35 d)
   which is a totally dry period, two, 7±2 days before parturition (−7 d)
   characterized by accelerated cell proliferation and pre-colostrum
   secretion (stage I of lactogenesis) and three, 3 days after parturition
   (+3 d) characterized by colostrum production (stage II of
   lactogenesis). Mammary biopsy was performed using Bard Magnum biopsy
   system (Bard Peripheral Vascular, Inc., Tempe, AZ, US). The collected
   samples were immediately frozen and stored in liquid nitrogen until
   further analysis.

RNA isolation

   To reduce the number of mammary samples needed for DGE analysis,
   samples from cows in the same group were pooled together for RNA
   isolation. Total RNA was extracted with Trizol reagent (Invitrogen,
   Carlsbad, CA, USA) in accordance with the manufacturer's instructions.
   The quantity and quality of the extracted RNA was determined using a
   spectrophotometer (Biophotometer Plus, Eppendorf, Germany), with RNA
   quality being evaluated by the absorbance ratio at 260 nm/280 nm.

High-throughput sequencing

   High-throughput sequencing of mRNA was performed on an Illumina HiSeq
   2000 platform. All sample preparation and sequencing procedures were
   executed according to the manufacturer's instructions (Illumina, Inc.,
   San Diego, CA, USA). Briefly, mRNA from 6 µg of total RNA extracted
   from samples at each time point was purified by adsorbing to Oligo (dT)
   magnetic beads, and then use Oligo (dT) as primer to synthesize the
   first and second-strand cDNA. The bead-bound cDNA was digested with
   restriction enzyme NlaIII, which recognized and cut off the CATG sites.
   The fragments apart from the 3′ cDNA fragments connected to Oligo (dT)
   beads were washed away and the Illumina adaptor 1 was ligated to the
   sticky 5′ end of the digested bead-bound cDNA fragments. The junction
   of Illumina adaptor 1 and CATG site was the recognition site of Mmel,
   which was a type of Endonuclease with separated recognition sites and
   digestion sites. It cut at 17 bp downstream of the CATG site, producing
   tags with adaptor 1. After removing 3′ fragments with magnetic beads
   precipitation, Illumina adaptor 2 was ligated to the 3′ ends of tags,
   acquiring tags with different adaptors of both ends to form a tag
   library. After 15 cycles of linear PCR amplification, 105 bp fragments
   are purified by 6% TBE PAGE Gel electrophoresis. After denaturation,
   the single-stranded cDNA was anchored on Illumina flowcells of the
   Illumina HiSeq 2000 system and sequenced. Raw data (tag sequences) were
   deposited in NCBI's Gene Expression Omnibus (GEO) database under
   submission number [48]GSE44796.

Tag mapping

   For the raw data, we filtered adaptor sequences, low quality tags (tags
   with unknown nucleotides N), empty reads and tags that were too short
   or too long, and tags with only one copy to get clean tags.

   A virtual library containing all the possible CATG+17 bases length
   sequences of the bovine gene sequences were constructed by BGI-shenzhen
   [49][17]–[50][19], using Bos taurus UMD 3.1 [51][20]. All clean tags
   were mapped to the reference sequences, and only upto 1-bp mismatch was
   considered. Clean tags that mapped to reference sequences from multiple
   genes were filtered out. The remaining clean tags were designated as
   unambiguous clean tags. The number of unambiguous clean tags for each
   gene was calculated and then normalized to TPM (number of transcripts
   per million clean tags).

Screening of differentially expressed genes

   A rigorous algorithm has been developed to identify genes expressed
   differentially between samples. Because every gene's expression
   occupies only a small part of the library, p(x) was within the Poisson
   distribution. The probability of a given gene being expressed equally
   between two samples can be calculated with the following equation
   [52][21]:
   graphic file with name pone.0070393.e001.jpg

   where the total number of clean tags in sample 1 is N[1], the total
   number of clean tags in sample 2 is N[2], and there are x tags in
   sample 1 and y tags in sample 2 for gene A. False Discovery Rate (FDR)
   analysis was applied to determine the threshold P value in multiple
   tests and analyses [53][22]. The significance of the difference in gene
   expression was judged using a threshold of FDR≤0.001 and fold-change
   value ≥2 (|log2 Ratio|≥1).

Sequencing data analysis

   The selected genes with significant modification in their expression
   were subjected to further analysis, which were GO enrichment analysis,
   pathway enrichment analysis and functional annotation clustering.

   GO enrichment analysis of functional significance applies
   hypergeometric test to map all differentially expressed genes (DEGs) to
   terms in Gene Ontology (GO) database (Version 1.1.1631) [54][23],
   looking for significantly enriched GO terms in DEGs comparing to the
   genome background. The calculating formula is [55][24]:
   graphic file with name pone.0070393.e002.jpg

   where N is the number of all genes with GO annotation; n is the number
   of DEGs in N; M is the number of all genes that are annotated to the
   certain GO terms; m is the number of DEGs in M. After Bonferroni
   correction, GO terms with corrected-p value ≤0.05 were significantly
   enriched in differentially expressed genes.

   Pathway enrichment analysis based on Kyoto Encyclopedia of Genes and
   Genomes (KEGG) database (Version 58) can identify significantly
   enriched biological pathways in DEGs compared with the whole genome
   background [56][25]. Enriched pathways were calculated using the same
   formula as that in GO analysis. Here N is the number of all genes with
   KEGG annotation, n is the number of DEGs in N, M is the number of all
   genes annotated in specific pathways, and m is the number of DEGs in M.
   Pathways with a Q value ≤0.05 were significantly enriched in
   differentially expressed genes.

   The differentially expressed genes detected in stage I and stage II of
   lactogenesis were further analyzed using DAVID functional annotation
   clustering tool [57][26]. This tool can provide a look at the internal
   relationships of the clustered terms.

Quantitative real-time PCR (qRT-PCR)

   QRT-PCR analysis was used to verify the DGE results. The RNA samples
   used for the qRT-PCR assay were both the same as for the DGE
   experiments and independent RNA extractions from biological replicates.
   13 candidate genes were selected and detected using qRT-PCR, including
   caseinβ, CIDEA, GLYCAM, BGN, ELL3, PAH, GP2, LOC525947, SLC14A1, SST,
   PNMT, caseinα-S2 and integrinβ. qRT-PCR was performed according to the
   TaKaRa manufacturer specifications (TaKaRa SYBR® PrimeScript™ RT-PCR
   Kit, Dalian, China). SYBR Green PCR cycling was denatured using a
   program of 95°C for 10 s, and 40 cycles of 95°C for 5 s and 60°C for 40
   s, and performed on an ABI 7500 instrument (Applied Biosystems, Foster
   City, CA, USA). Data were reported as values normalized to the
   housekeeping gene β-actin, and they were subjected to one-way ANOVA
   analysis using Statistical Analysis Systems statistical software
   package (Version 8e, SAS Institute, Cary, NC, USA). Means were
   considered significantly different at p<0.05.

Results

Analysis of the three cDNA libraries

   Characterization of the three cDNA libraries (−35 d, −7 d and +3 d) was
   listed in [58]Table 1. Sequencing depths of 6,203,209, 5,857,755 and
   6,115,741 tags were identified in the three libraries, which included
   192,928, 170,427 and 168,233 distinct tags, respectively. As the copy
   number of the tags increased, the number of distinct tags decreased.
   Saturation analysis was performed to check whether the number of
   detected genes continues to increase as the number of tags sequenced
   (total tag number) increases. The results were shown in [59]Figure 1.
   When the number of sequenced tags reached 2M or more, the number of
   detected genes almost ceased increasing. There were 6.2M, 5.8M and 6.1M
   tags in the three libraries; thus, all three were sequenced to
   saturation, producing a full representation of the transcripts present
   in each of the three stages.

Table 1. Distribution of the tags sequenced from the three libraries.

   Tag library
   −35 d −7 d +3 d
   Total tag number 6,203,209 5,857,755 6,115,741
   Total clean tag number 192,928 170,427 168,233
   Distinct clean tag number 476,719 429,231 430,516
   Tag copy number <2 273,198 250,858 254,875
   Tag copy number ≥2 (clean tags) 192,928 (100%) 170,427 (100%) 168,233
   (100%)
   Tag copy number >5 79,694 (41.31%) 67,230 (39.45%) 63,337 (37.65%)
   Tag copy number >10 50,966 (26.42%) 42,877 (25.16%) 39,889 (23.71%)
   Tag copy number >20 31,557 (16.36%) 26,947 (15.81%) 24,676 (14.67%)
   Tag copy number >50 15,682 (8.13%) 14,034 (8.23%) 12,451 (7.40%)
   Tag copy number >100 8,427 (4.37%) 7,864 (4.61%) 6,693 (3.98%)
   [60]Open in a new tab

Figure 1. Saturation analysis of −35 d (A), −7 d (B) and +3 d (C) libraries
sequencing.

   Figure 1
   [61]Open in a new tab

   “—” represents genes mapped by all clean tags. “---” represents genes
   mapped by unambiguous clean tags. The number of detected genes
   continued increasing as the total tag number increased. When the number
   of total tags reached 2M or more, the number of detected genes almost
   ceased increasing.

Analysis of tag annotation

   A reference database of the bovine genome includes 35,945 transcripts
   [62][20]. A total of 193,547 reference tags with 94,396 (48.77%)
   unambiguous reference tags were obtained. All clean tags in the three
   cDNA libraries were mapped to the reference tags, and the results were
   shown in [63]Figure 2. There were 59,537 (30.86%), 48,944 (28.72%) and
   48,742 (28.98%) distinct tags matched to the reference genes in the
   three libraries (−35 d, −7 d and +3 d), respectively. Among these tags,
   there were 8.99%, 9.60% and 9.33% tags, respectively, matched to the
   anti-sense strand of the genes, suggesting that the antisense strand of
   these genes also had transcripts and that these genes might therefore
   have sense-antisense regulation. The unmatched tags were then mapped to
   the bovine genome, and 37.90%, 40.08% and 39.40% of tags, respectively,
   were matched to the genomic sequences in the three libraries. These
   tags might represent non-annotated genes or could possibly be derived
   from intergenic regions not encoding any transcripts.

Figure 2. Summary of the distinct tag-to-gene mapping data in the three
libraries.

   Figure 2
   [64]Open in a new tab

Detection of differentially expressed genes and GO enrichment analysis

   At the cutoff criteria of FDR≤0.001 and |log2 Ratio|≥1, many genes were
   differentially expressed; these results were shown in [65]Figure 3. The
   details of these differentially expressed genes were supplied in Excel
   S1, S2 and S3. A total of 812 genes were significantly differentially
   expressed at −7 d compared with −35 d (stage I), accounting for 9.70%
   of transcripts in the −7 d cDNA library. There were 234 (28.80%) genes
   upregulated (in red) and 578 (71.20%) genes downregulated (in green).
   The number of genes having greater than a five-fold difference
   accounted for 0.45% of the total differentially expressed genes
   ([66]Figure 4). A total of 1,189 (14.20%) genes were significantly
   differentially expressed, with 274 (23.00%) genes upregulated and 915
   (77.00%) genes downregulated at +3 d compared with −7 d (stage II). The
   number of genes having greater than a five-fold difference accounted
   for 0.41% of the total differentially expressed genes. Moreover, there
   were 1,672 genes significantly differentially expressed at +3 d
   compared with −35 d. There were 209 (12.5%) upregulated genes, and
   1,463 (87.5%) downregulated genes. Among those genes, there were 11
   (0.65%) genes having greater than a five-fold difference.

Figure 3. Pairwise comparisons of the level of gene expression between
libraries.

   Figure 3
   [67]Open in a new tab

   To compare the level of gene expression between two libraries, each
   library was normalized to TPM. (A) Differentially expressed genes at −7
   d compared with −35 d. (B) Differentially expressed genes at +3 d
   compared with −7 d. (C) Differentially expressed genes at +3 d compared
   with −35 d. Red dots represent upregulated transcripts. Green dots
   represent downregulated transcripts. Blue dots represent transcripts
   that did not change significantly. The parameters FDR≤0.001 and |log2
   Ratio|≥1 were used as the threshold to judge the significance of the
   difference in gene expression. See supplementary data 1, 2 and 3 for
   details of these differentially expressed genes.

Figure 4. Differentially expressed tags at the onset of lactation.

   Figure 4
   [68]Open in a new tab

   (A) Differentially expressed tags at −7 d compared with −35 d. (B)
   Differentially expressed tags at +3 d compared with −7 d. (C)
   Differentially expressed tags at +3 d compared with −35 d. The “x”
   represents the fold-change of differentially expressed unique tags. The
   “y” axis represents the number of unique tags (log10). The red region
   represents differentially accumulating unique tags with a 5-fold
   difference between libraries. The green and blue regions represent
   unique tags that are upregulated and downregulated for more than 5
   fold.

   The differentially expressed genes were analyzed in the context of gene
   ontology (GO) biological processes. This analysis revealed that the
   upregulated genes in DEGs were mainly associated with the cell cycle
   and lipid metabolism, and the downregulated genes were mainly
   associated with the immune response and biological adhesion at −7 d
   compared with −35 d. And the upregulated genes in DEGs at +3 d compared
   with −7 d were mainly associated with the immune response, and the
   downregulated genes were mainly associated with the cell cycle phase.
   The most significantly enriched biological processes in DEGs were
   metabolic processes at +3 d compared with −35 d.

   Moreover, during the onset of lactation, some genes were completely
   turned on or off. There were 305 genes turned on at −7 d compared with
   −35 d. Gene ontology analysis showed that these genes were associated
   with signal transduction (GABRB1; GPRC5A; GPR107; RAB9B), transport
   (SLC5A1; SLC13A3; KCNC1; KCNAB3; SLC38A5), the cell cycle (CDCA5;
   HORMAD1; FSD1; MPHOSPH9), and other processes. In addition, 12 genes
   encoding hypothetical proteins (LOC100335323, LOC100336058,
   LOC100335520) were also identified in stage I of lactogenesis. 604
   genes were turned off at −7 d compared with −35 d. Gene ontology
   analysis showed that these genes were associated with ion transport
   (SLC26A2; SCN4B; TMEM163; SLC5A12; SLC13A4; SLC39A11; SLCO2A1) and
   protein and lipid metabolism (PREP; caspase-15; LPGAT1; PLCB4; EPT1).
   There were 421 genes turned on at +3 d compared with −7 d, including
   genes involved in transport (SLC5A2; SLC26A2; KCNH2), the immune
   response (MICB; CCL28; IL8), and signal transduction (PIK3R3; AKAP4).
   425 genes were turned off at +3 d compared with −7 d, including cell
   cycle-related genes (CDC25C; SKA1; FSD1) and signal transduction genes
   (GABRB1; GPRC5A; GPR37).

Pathway enrichment analysis of differentially expressed genes

   After pathway enrichment analysis, the top 10 significantly enriched
   pathways were listed in [69]Tables 2, [70]3 and [71]4. The results
   suggested that many of the pathways enriched in DEGs were associated
   with the immune response at −7 d compared with −35 d, while many of the
   pathways enriched in DEGs were associated with the cell cycle and
   metabolism at +3 d compared with −7 d. And the most enriched pathway in
   DEGs was metabolic pathways at +3 d compared with −35 d.

Table 2. Top 10 significantly enriched pathways at −7 d compared with −35 d.

   Pathway DEGs with pathway annotation Q value Pathway ID
   Phagosome 41 (6.23%) 2.87e-10 ko04145
   Complement and coagulation cascades 22 (3.34%) 3.12e-04 ko04610
   Lysosome 23 (3.5%) 3.12e-04 ko04142
   Rheumatoid arthritis 18 (2.74%) 5.85e-04 ko05323
   Hematopoietic cell lineage 17 (2.58%) 1.60e-03 ko04640
   Staphylococcus aureus infection 19 (2.89%) 1.71e-03 ko05150
   Cell adhesion molecules (CAMs) 26 (3.95%) 1.84e-03 ko04514
   Metabolic pathways 90 (13.68%) 4.16e-03 ko01100
   PPAR signaling pathway 15 (2.28%) 9.70e-03 ko03320
   Bladder cancer 9 (1.37%) 9.70e-03 ko05219
   [72]Open in a new tab

Table 3. Top 10 significantly enriched pathways at +3 d compared with −7 d.

   Pathway                   DEGs with pathway annotation Q value  Pathway ID
   Metabolic pathways        172 (18.03%)                 3.74e-15 ko01100
   Proteasome                18 (1.89%)                   6.43e-10 ko03050
   DNA replication           15 (1.57%)                   5.50e-08 ko03030
   Base excision repair      14 (1.47%)                   3.04e-06 ko03410
   Cell cycle                30 (3.14%)                   6.36e-06 ko04110
   Pyrimidine metabolism     23 (2.41%)                   4.32e-05 ko00240
   Purine metabolism         33 (3.46%)                   5.13e-05 ko00230
   Oxidative phosphorylation 25 (2.62%)                   5.90e-05 ko00190
   Propanoate metabolism     12 (1.26%)                   9.46e-05 ko00640
   Mismatch repair           9 (0.94%)                    2.17e-04 ko03430
   [73]Open in a new tab

Table 4. Top 10 significantly enriched pathways in differentially expressed
genes at +3 d compared with −35 d.

   Pathway                    DEGs with pathway annotation Q value  Pathway ID
   Metabolic pathways         205 (15.67%)                 1.32e-11 ko01100
   Oxidative phosphorylation  38 (2.91%)                   2.35e-08 ko00190
   DNA replication            16 (1.22%)                   4.23e-07 ko03030
   Proteasome                 17 (1.3%)                    4.23e-07 ko03050
   Base excision repair       16 (1.22%)                   2.99e-06 ko03410
   Cell cycle                 35 (2.68%)                   2.16e-05 ko04110
   Phagosome                  48 (3.67%)                   2.16e-05 ko04145
   Parkinson's disease        35 (2.68%)                   2.95e-05 ko05012
   Nucleotide excision repair 16 (1.22%)                   3.78e-04 ko03420
   Huntington's disease       45 (3.44%)                   5.80e-04 ko05016
   [74]Open in a new tab

Analysis of milk protein genes and lipogenic genes

   The expression level of many milk protein genes increased significantly
   during stage I of lactogenesis and kept increasing during stage II.
   These genes include casein kappa (CSN3), lactalbumin alpha (LALBA),
   casein beta (CSN2), casein alpha s1 (CSN1S1), and casein alpha-S2
   (CSN1S2). The fold changes of these genes were listed in [75]Table 5.
   Furthermore, many lipogenic genes changed significantly at the onset of
   lactation. These genes were associated with mammary fatty acids uptake
   from the blood (LPL, CD36), intracellular fatty acids trafficking
   (FABP3), long-chain (ACSL1) and short-chain (ACSS2) intracellular fatty
   acids activation, de novo fatty acids synthesis (ACACA, FASN),
   desaturation (SCD), triacylglycerol synthesis (AGPAT6, GPAM), lipid
   droplet formation (BTN1A1, XDH), ketone body utilization (BDH1), and
   transcription regulation (INSIG1, PPARGC1A). The fold changes of these
   genes were listed in [76]Table 6.

Table 5. Expression changes of milk protein genes during the two stages of
lactogenesis.

   Gene Description TPM (−35 d) TPM (−7 d) TPM (3 d) log[2] Ratio (−7/−35)
   log[2] Ratio (+3/−7)
   gi|31341749|ref|[77]NM_174528.2| Casein alpha –S2 (CSN1S2) 22.49 304.12
   10389.3 3.76 5.09
   gi|31342165|ref|[78]NM_174378.2| Lactalbumin alpha (LALBA) 3.72 198.63
   5224.76 5.74 4.72
   gi|31341348|ref|[79]NM_181029.2| Casein alpha s1 (CSN1S1) 578.92
   10599.41 80620.84 4.19 2.93
   gi|31341343|ref|NM_181008.2| Casein beta (CSN2) 735.32 12496.16
   64462.97 4.09 2.37
   gi|27881411|ref|NM_174294.1| Casein kappa (CSN3) 1390.84 4109.25
   13441.57 1.56 1.71
   [80]Open in a new tab

   Note: TPM (Transcripts Per Million clean tags) is a standardized
   indicator that specifies the number of transcript copies in every one
   million clean tags.

Table 6. Expression changes of lipogenic genes during the two stages of
lactogenesis.

   Gene Description TPM (−35 d) TPM (−7 d) TPM (+3 d) log[2] Ratio
   (−7/−35) log[2] Ratio (+3/−7)
   gi|115497163|ref|[81]NM_001075120.1| lipoprotein lipase (LPL) 35.17
   38.26 172.26 0.12 2.17[82]^*
   gi|31343049|ref|NM_174010.2| CD36 molecule (thrombospondin receptor)
   (CD36) 156.9 77.24 272.33 −1.02[83]^* 1.82[84]^*
   gi|31342355|ref|[85]NM_174313.2| fatty acid binding protein 3, muscle
   and heart (mammary-derived growth inhibitor) (FABP3) 1.52 46.31 156.52
   4.93[86]^* 1.76[87]^*
   gi|155372064|ref|[88]NM_001101169.1| solute carrier family 27 (fatty
   acid transporter) member 6 (SLC27A6) 9.64 36.12 50.98 1.91[89]^* 0.50
   gi|115497269|ref|[90]NM_001076085.1| acyl-CoA synthetase long-chain
   family member 1 (ACSL1) 96.37 174.68 115.64 0.86 −0.60
   gi|157427805|ref|[91]NM_001105339.1| acyl-CoA synthetase short-chain
   family member 2 (ACSS2) 30.27 78.85 32.33 1.38[92]^* −1.29[93]^*
   gi|31342550|ref|[94]NM_174224.2| acetyl-CoA carboxylase alpha (ACACA)
   1.86 2.5 2.74 0.43 0.13
   gi|60592789|ref|NM_001012669.1| fatty acid synthase (FASN) 660.42
   2804.46 652.95 2.09[95]^* −2.10[96]^*
   gi|148540093|ref|[97]NM_173959.4| stearoyl-CoA desaturase
   (delta-9-desaturase) (SCD) 74.73 796.32 471.62 3.41[98]^* −0.76
   gi|139948314|ref|[99]NM_001083669.1| 1-acylglycerol-3-phosphate
   O-acyltransferase 6 (lysophosphatidic acid acyltransferase, zeta)
   (AGPAT6) 41.59 32.72 44.48 −0.35 0.44
   gi|59676567|ref|[100]NM_001012282.1| glycerol-3-phosphate
   acyltransferase, mitochondrial (GPAM) 77.27 113.53 107.77 0.56 −0.075
   gi|148232524|ref|[101]NM_205793.2| diacylglycerol O-acyltransferase 2
   (DGAT2) 12.17 29.5 2.57 1.28[102]^* −3.52[103]^*
   gi|31341801|ref|[104]NM_174508.2| butyrophilin subfamily 1 member A1
   (BTN1A1) 4.4 42.55 85.87 3.27[105]^* 1.01[106]^*
   gi|31343144|ref|[107]NM_173972.2| xanthine dehydrogenase (XDH) 299.44
   827.07 1390.57 1.47[108]^* 0.75
   gi|77736146|ref|NM_001034600.1| 3-hydroxybutyrate dehydrogenase, type 1
   (BDH1) 0.51 2.5 2.91 2.29 0.22
   gi|118150923|ref|[109]NM_001077909.1| insulin induced gene 1 (INSIG1)
   40.24 94.94 65 1.24[110]^* −0.55
   gi|164519007|ref|[111]NM_001113302.1| sterol regulatory element binding
   transcription factor 1 (SREBF1) 57.99 147.32 92.03 1.35[112]^* −0.68
   [113]Open in a new tab

   Note:
   ^*

   represents significant difference at the cutoff criteria of FDR≤0.001
   and |log2 Ratio|≥1.

Analysis of differentially expressed genes throughout the onset of lactation

   There were 231 genes differentially expressed in the two stages of
   lactogenesis. Among these genes, 92 genes were downregulated in stage I
   and upregulated in stage II, 111 genes were upregulated in stage I and
   downregulated in stage II, and only 11 genes continued to increase in
   both stages. After analysis of functional annotation clustering, 92
   genes fell into 24 clusters. Each term in Cluster 4 was associated with
   the immune response, and each in cluster two was associated with the
   extracellular matrix. The top 5 enriched clusters were listed in
   [114]Table 7. After pathway enrichment analysis, the ECM-Receptor
   interaction was enriched ([115]Figure 5). In addition, functional
   annotation clustering of the 111 genes that were upregulated initially
   and then downregulated fell into 19 clusters; the top 5 enriched
   clusters were listed in [116]Table 8. We can find that all terms in
   Cluster 1 were associated with cell cycle.

Table 7. Functional annotation clustering of genes downregulated in stage I
and then upregulated in stage II of lactogenesis.

   Category Terms N Percent(%) P-Value
   Cluster 1 Enrichment Score: 12.80
   GOTERM_CC_FAT GO:0005576∼extracellular region 36 41.38 2.88E-17
   UP_SEQ_FEATURE signal peptide 37 42.53 6.67E-15
   SP_PIR_KEYWORDS signal 37 42.53 1.07E-14
   SP_PIR_KEYWORDS Secreted 29 33.33 9.55E-14
   SP_PIR_KEYWORDS disulfide bond 32 36.78 9.52E-13
   UP_SEQ_FEATURE glycosylation site:N-linked (GlcNAc…) 34 39.08 3.16E-12
   SP_PIR_KEYWORDS glycoprotein 35 40.23 6.54E-12
   UP_SEQ_FEATURE disulfide bond 29 33.33 9.64E-11
   Cluster 2 Enrichment Score: 5.50
   GOTERM_CC_FAT GO:0044421 extracellular region part 20 22.99 8.79E-10
   GOTERM_CC_FAT GO:0031012 extracellular matrix 11 12.64 1.96E-06
   SP_PIR_KEYWORDS extracellular matrix 7 8.05 1.36E-04
   GOTERM_CC_FAT GO:0005578 proteinaceous extracellular matrix 8 9.20
   4.14E-04
   Cluster 3 Enrichment Score: 5.46
   INTERPRO IPR014716:Fibrinogen, alpha/beta/gamma chain, C-terminal
   globular, subdomain 1 6 6.90 2.87E-07
   INTERPRO IPR002181:Fibrinogen, alpha/beta/gamma chain, C-terminal
   globular 6 6.90 7.05E-07
   SMART SM00186:FBG 6 6.90 3.55E-06
   UP_SEQ_FEATURE domain:Fibrinogen C-terminal 4 4.60 2.02E-04
   Cluster 4 Enrichment Score: 3.68
   GOTERM_BP_FAT GO:0002526 acute inflammatory response 7 8.05 5.29E-07
   GOTERM_BP_FAT GO:0006954 inflammatory response 8 9.20 8.54E-06
   GOTERM_BP_FAT GO:0002541 activation of plasma proteins involved in
   acute inflammatory response 5 5.75 8.59E-06
   GOTERM_BP_FAT GO:0006956 complement activation 5 5.75 8.59E-06
   GOTERM_BP_FAT GO:0006957 complement activation, alternative pathway 4
   4.60 1.21E-05
   SP_PIR_KEYWORDS complement alternate pathway 4 4.60 2.82E-05
   GOTERM_BP_FAT GO:0006959 humoral immune response 5 5.75 3.19E-05
   GOTERM_BP_FAT GO:0009611 response to wounding 9 10.34 3.22E-05
   GOTERM_BP_FAT GO:0045087 innate immune response 6 6.90 7.05E-05
   GOTERM_BP_FAT GO:0051605 protein maturation by peptide bond cleavage 5
   5.75 1.48E-04
   GOTERM_BP_FAT GO:0050778 positive regulation of immune response 6 6.90
   2.14E-04
   GOTERM_BP_FAT GO:0002253 activation of immune response 5 5.75 3.40E-04
   SP_PIR_KEYWORDS innate immunity 5 5.75 3.93E-04
   GOTERM_BP_FAT GO:0016485 protein processing 5 5.75 5.02E-04
   KEGG_PATHWAY bta04610:Complement and coagulation cascades 5 5.75
   5.36E-04
   GOTERM_BP_FAT GO:0051604 protein maturation 5 5.75 6.23E-04
   GOTERM_BP_FAT GO:0006952 defense response 8 9.20 0.00
   GOTERM_BP_FAT GO:0048584 positive regulation of response to stimulus 6
   6.90 0.00
   GOTERM_BP_FAT GO:0002252 immune effector process 5 5.75 0.00
   GOTERM_BP_FAT GO:0002684 positive regulation of immune system process 6
   6.90 0.00
   SP_PIR_KEYWORDS immune response 5 5.75 0.00
   GOTERM_BP_FAT GO:0006955 immune response 7 8.05 0.02
   GOTERM_BP_FAT GO:0006508 proteolysis 9 10.34 0.07
   Cluster 5 Enrichment Score: 2.84
   KEGG_PATHWAY bta04512:ECM-receptor interaction 7 8.05 2.30E-06
   INTERPRO IPR013032:EGF-like region, conserved site 9 10.34 7.39E-06
   INTERPRO IPR006210:EGF-like 7 8.05 1.20E-04
   SMART SM00181:EGF 7 8.05 6.99E-04
   INTERPRO IPR013111:EGF, extracellular 4 4.60 0.00
   GOTERM_BP_FAT GO:0007155 cell adhesion 9 10.34 0.00
   GOTERM_BP_FAT GO:0022610 biological adhesion 9 10.34 0.00
   KEGG_PATHWAY bta04510:Focal adhesion 6 6.90 0.00
   SP_PIR_KEYWORDS cell adhesion 6 6.90 0.00
   GOTERM_BP_FAT GO:0007160 cell-matrix adhesion 3 3.45 0.03
   GOTERM_BP_FAT GO:0031589 cell-substrate adhesion 3 3.45 0.03
   INTERPRO IPR002035:von Willebrand factor, type A 3 3.45 0.04
   SMART SM00327:VWA 3 3.45 0.08
   [117]Open in a new tab

Figure 5. The ECM-Receptor interaction pathway is transcriptionally regulated
throughout the onset of lactation.

   [118]Figure 5
   [119]Open in a new tab

   Red nodes represent genes downregulated during stage I of lactogenesis
   and then upregulated during stage II of lactogenesis.

Table 8. Functional annotation clustering of genes upregulated during stage I
and then downregulated during stage II of lactogenesis.

   Category Term N Percent(%) P-Value
   Cluster 1 Enrichment Score: 2.48
   GOTERM_BP_FAT GO:0007049 cell cycle 11 10.68 7.85E-06
   GOTERM_BP_FAT GO:0022402 cell cycle process 9 8.74 2.51E-05
   GOTERM_BP_FAT GO:0022403 cell cycle phase 8 7.77 3.66E-05
   GOTERM_BP_FAT GO:0000279 M phase 7 6.80 9.29E-05
   GOTERM_BP_FAT GO:0000278 mitotic cell cycle 7 6.80 1.03E-04
   GOTERM_BP_FAT GO:0000280 nuclear division 6 5.83 1.77E-04
   GOTERM_BP_FAT GO:0007067 mitosis 6 5.83 1.77E-04
   SP_PIR_KEYWORDS cell division 7 6.80 2.23E-04
   GOTERM_BP_FAT GO:0000087 M phase of mitotic cell cycle 6 5.83 2.24E-04
   GOTERM_BP_FAT GO:0048285 organelle fission 6 5.83 2.37E-04
   GOTERM_BP_FAT GO:0051301 cell division 6 5.83 0.00
   GOTERM_CC_FAT GO:0000775 chromosome, centromeric region 5 4.85 0.00
   SP_PIR_KEYWORDS cell cycle 7 6.80 0.00
   SP_PIR_KEYWORDS mitosis 5 4.85 0.00
   GOTERM_CC_FAT GO:0005694 chromosome 8 7.77 0.00
   GOTERM_CC_FAT GO:0044427 chromosomal part 7 6.80 0.01
   GOTERM_CC_FAT GO:0005819 spindle 4 3.88 0.01
   GOTERM_CC_FAT GO:0015630 microtubule cytoskeleton 6 5.83 0.03
   GOTERM_BP_FAT GO:0007017 microtubule-based process 4 3.88 0.06
   GOTERM_CC_FAT GO:0000793 condensed chromosome 3 2.91 0.06
   SP_PIR_KEYWORDS cytoskeleton 5 4.85 0.10
   GOTERM_CC_FAT GO:0043228 non-membrane-bounded organelle 13 12.62 0.18
   GOTERM_CC_FAT GO:0043232 intracellular non-membrane-bounded organelle
   13 12.62 0.18
   GOTERM_CC_FAT GO:0005874 microtubule 3 2.91 0.19
   GOTERM_CC_FAT GO:0044430 cytoskeletal part 6 5.83 0.21
   GOTERM_CC_FAT GO:0005856 cytoskeleton 7 6.80 0.30
   SP_PIR_KEYWORDS cytoplasm 10 9.71 0.69
   Cluster 2 Enrichment Score: 2.34
   GOTERM_MF_FAT GO:0048037∼cofactor binding 9 8.74 3.42E-05
   GOTERM_MF_FAT GO:0050662∼coenzyme binding 5 4.85 0.01
   GOTERM_BP_FAT GO:0055114∼oxidation reduction 6 5.83 0.261
   Cluster 3 Enrichment Score: 1.71
   GOTERM_MF_FAT GO:0030554∼adenyl nucleotide binding 20 19.42 4.08E-05
   GOTERM_MF_FAT GO:0001883∼purine nucleoside binding 20 19.42 4.63E-05
   GOTERM_MF_FAT GO:0001882∼nucleoside binding 20 19.42 5.00E-05
   GOTERM_MF_FAT GO:0032559∼adenyl ribonucleotide binding 17 16.50
   7.63E-04
   GOTERM_MF_FAT GO:0017076∼purine nucleotide binding 20 19.42 8.58E-04
   GOTERM_MF_FAT GO:0005524∼ATP binding 16 15.53 0.00
   SP_PIR_KEYWORDS atp-binding 13 12.62 0.00
   GOTERM_MF_FAT GO:0000166∼nucleotide binding 21 20.39 0.00
   GOTERM_MF_FAT GO:0032553∼ribonucleotide binding 17 16.50 0.01
   GOTERM_MF_FAT GO:0032555∼purine ribonucleotide binding 17 16.50 0.01
   INTERPRO IPR017441:Protein kinase, ATP binding site 6 5.83 0.03
   INTERPRO IPR000719:Protein kinase, core 6 5.83 0.06
   SP_PIR_KEYWORDS nucleotide-binding 12 11.65 0.06
   GOTERM_MF_FAT GO:0004674∼protein serine/threonine kinase activity 5
   4.85 0.11
   GOTERM_MF_FAT GO:0004672∼protein kinase activity 6 5.83 0.15
   SP_PIR_KEYWORDS serine/threonine-protein kinase 4 3.88 0.15
   GOTERM_BP_FAT GO:0006468∼protein amino acid phosphorylation 6 5.83 0.18
   INTERPRO IPR008271:Serine/threonine protein kinase, active site 4 3.88
   0.19
   SMART SM00220:S_TKc 3 2.91 0.22
   INTERPRO IPR017442:Serine/threonine protein kinase-related 4 3.88 0.24
   GOTERM_BP_FAT GO:0016310∼phosphorylation 6 5.83 0.27
   INTERPRO IPR002290:Serine/threonine protein kinase 3 2.91 0.29
   SP_PIR_KEYWORDS kinase 4 3.88 0.42
   GOTERM_BP_FAT GO:0006793∼phosphorus metabolic process 6 5.83 0.44
   GOTERM_BP_FAT GO:0006796∼phosphate metabolic process 6 5.83 0.44
   Cluster 4 Enrichment Score: 1.37
   GOTERM_BP_FAT GO:0006259∼DNA metabolic process 7 6.80 0.01
   GOTERM_BP_FAT GO:0006281∼DNA repair 4 3.88 0.06
   GOTERM_BP_FAT GO:0033554∼cellular response to stress 5 4.85 0.07
   GOTERM_BP_FAT GO:0006974∼response to DNA damage stimulus 4 3.88 0.10
   Cluster 5 Enrichment Score: 1.28
   GOTERM_BP_FAT GO:0008203∼cholesterol metabolic process 3 2.91 0.03
   GOTERM_BP_FAT GO:0016125∼sterol metabolic process 3 2.91 0.04
   GOTERM_BP_FAT GO:0008202∼steroid metabolic process 3 2.91 0.11
   [120]Open in a new tab

Quantitative real-time PCR (qRT-PCR) confirmation

   The results of qRT-PCR were shown in [121]Figure 6, the expression
   patterns of 13 genes showed a general agreement with the results of the
   DGE experiments, suggesting that DGE was an efficient and accurate
   strategy for the detection of differentially expressed genes.

Figure 6. Comparison of the results of 13 genes obtained using DGE and
qRT-PCR methodology.

   Figure 6
   [122]Open in a new tab

Discussion

   The major goal of this study was to explore genome-wide gene expression
   profiles during the onset of lactation and to provide information for
   understanding the underlying molecular mechanisms. Three libraries were
   constructed from mammary tissues at three different stages (−35 d, −7 d
   and +3 d relative to parturition). Saturation analysis of the sequences
   indicated that the three libraries had reached saturation and were
   therefore complete assessments of all transcripts present in the
   libraries. The three libraries contained 192,928, 170,427 and 168,233
   distinct clean tags, respectively. However, a reference database
   includes 35,945 transcripts. Theoretically, a tag should be generated
   by NlaIII digestion from the 3′-most end of a transcript. However, many
   transcripts contain more than one CATG site, so tags from other NlaIII
   sites were also generated in our libraries. Because only one tag can be
   generated per transcript from a given NlaIII site in a cDNA, these
   extra NlaIII tags represent genes redundantly present in the expression
   profile. Thus, the number of unique tags generated was greater than
   that of the annotated bovine genome.

   To validate the DGE method, the levels of 13 genes were analyzed using
   qRT-PCR. Although the differences in the expression of some genes did
   not match the magnitude of those detected using the DGE method, the
   trends of upregulation and downregulation were similar. Compared with
   microarray and qRT-PCR methods, the sequencing method has been
   documented to be more sensitive for the estimation of gene expression,
   especially for low-abundance transcripts [123][17].

Changes in mammary epithelial cells and the extracellular matrix at the onset
of lactation

   During the transition from pregnancy to lactation in the dairy cow, the
   mammary gland undergoes dramatic functional and metabolic changes,
   including the morphogenesis of mammary ducts during early pregnancy and
   differentiation of the mammary alveolus during late pregnancy
   [124][27], [125][28]. This study had found that many genes associated
   with the cell cycle were upregulated in stage I and downregulated in
   stage II, including cell division cycle associated 3 (CDCA3), Protein
   FAM83D (FAM83D), coiled-coil domain containing 99 (CCDC99), and others.
   These results suggested that cell proliferation occurred in late
   pregnancy and almost ceased after parturition, which was consistent
   with the early study [126][13]. Genes encoding extracellular matrix
   proteins were downregulated in stage I and upregulated in stage II. In
   addition, the ECM-Receptor interaction pathway was significantly
   enriched, and the differentially expressed genes in this pathway were
   downregulated and then later upregulated. These results demonstrated
   that, in stage I, the communication between cells and the extracellular
   matrix became weak with the proliferating of mammary epithelial cells
   and decreasing of the extracellular matrix, which recovered in stage
   II. A study in mice found that the communication between mammary
   epithelial cells and their environment became weak during the lactation
   period [127][11].

Changes in genes associated with the immune response

   At −7 d compared with −35 d, genes associated with the immune response
   accounted for 14.5% of differentially expressed genes, and genes
   associated with the defense response accounted for 8.20% of
   differentially expressed genes. Most of these genes were downregulated;
   for example, Chemokine (C-C motif) ligand 28 (CCL28) was completed
   turned off at −7 d. CCL28 is selectively expressed in certain mucosal
   tissues such as the exocrine glands, trachea, and colon and has a
   potent antimicrobial activity against Candida albicans, Gram-negative
   bacteria, and Gram-positive bacteria [128][29], [129][30]. Udder
   infections and mastitis are major problems for the dairy industry
   throughout the world, and the yearly costs are substantial. The risk of
   udder infection is the highest during the drying-off period and around
   the time of parturition [130][31], [131][32]. Adequate immune function
   is essential for the defense against udder infections. These genes
   associated with the immune response might be important for the immune
   function of the mammary gland and thus could be candidate genes for
   improving immune function.

Expression profiles and regulation of milk protein genes

   In this study, many milk protein genes were upregulated in stage I of
   lactogenesis. It has been shown previously in mouse [132][33], rat
   [133][34], and rabbit [134][35] that the expression of milk protein
   genes starts in early to mid-pregnancy, increases throughout pregnancy
   and reaches a plateau in late pregnancy and early lactation. However,
   Finucane et al [135][12] found that milk protein genes did not show
   significant changes in expression from late pregnancy to early
   lactation but that protein expression of glucose transporter GLUT1 did
   increase. The regulation of protein synthesis, particularly
   translation, in all mammalian tissues appears to be under control of
   the mTOR pathway [136][9]. Recently, studies in rodents and ruminants
   found that the mTOR pathway was very important for milk protein
   synthesis [137][36]–[138][40]. However, our results indicated that AKT
   (−1.10*) and PI3K (−1.20*) were downregulated significantly at −7 d
   compared with −35 d, which were upstream regulators of mTOR. In
   addition, some amino acid transporters were also downregulated during
   stage I, e.g. SLC43A2 (−1.85*), SLC3A2 (−1.58*) and SLC14A1 (−1.04*).
   In the mouse oocyte, mRNAs can be stored and then activated at the
   proper time [139][41]. Further study is required to investigate whether
   mRNA storage occurs in the bovine mammary gland.

   The main regulator of milk protein expression in non-ruminant mammary
   glands appears to be the Jak-Stat5 signaling pathway [140][42]. In
   addition to protein synthesis, STAT5 is important for mammary gland
   development [141][43]. In bovines, STAT5 responds to prolactin and
   other lactogenic growth factors, and its activity increases during
   lactation [142][44], [143][45]. However, when compared with the rodent
   mammary gland, the role of bovine STAT5 in controlling milk protein
   expression through the Jak-Stat5 signaling pathway appears to be
   considerably weaker [144][46]. Our data appear to support a minor role
   of Jak-Stat5 signaling in milk protein synthesis, implied by the lack
   of any changes in the expression of PRLR and STAT5B during lactation.
   However, STAT5B activity is mostly regulated by its phosphorylation
   status, and this appears to be already regulated at the onset of
   lactation. Further study is required to investigate whether STAT5
   phosphorylation increases at the onset of lactation.

Lipid metabolism and its regulation at the onset of lactation

   It has been speculated that there are two modes of regulation to
   control fatty acid synthesis in the mammary gland of the lactating
   mouse: the well-known SREBF1 system and a novel mechanism that acts at
   the posttranscriptional level in SCAP deletion mice fed a high-fat
   diet, leading to alterations to enzyme protein [145][47]. In this
   study, sterol regulatory element binding transcription factor 1
   (SREBF1) and cofactor insulin induced gene 1 (INSIG1) increased
   2.5-fold and 2.4-fold at −7 d compared with −35 d, respectively. It is
   possible that SREBF1 and cofactor INSIG1 regulate fat synthesis at the
   onset of lactation [146][48]. It has been postulated that, with the
   activation of secretion, SREBP-1 and its congener SREBP-2, the
   activators of cholesterol biosynthesis, are shuttled from the
   endoplasmic reticulum to the Golgi apparatus. Once in the Golgi, they
   are activated by proteolytic cleavage of a cytoplasmic fragment, a bHLH
   transcription factor that travels to the nucleus, directly activating
   the genes required for the synthesis of fatty acids and cholesterol
   [147][49]. This mechanism of action has been well studied in liver
   [148][50] and adipose tissues [149][51], and the promoters of the
   activated genes include SREs (sterol response elements) as well as
   binding sites for NF-Y, USF, SP1, and SP3 [150][52]–[151][54].

   It is well known that the onset of lactation includes two different
   stages. And in stage I of lactogenesis, the mammary epithelial cells
   proliferate and differentiate, the extracellular matrix diminishes, and
   the communication between the two becomes weak. As milk protein genes
   and lipogenesis genes are upregulated, genes associated with the immune
   response are downregulated. In stage II of lactogenesis, the mammary
   gland restores the communication between mammary epithelial cells and
   the extracellular matrix. The main result of this process is the
   synthesis of milk using nutrients from the blood.

Supporting Information

   File S1

   The details of differentially expressed genes at −7 d compared with −35
   d.

   (XLS)
   [152]Click here for additional data file.^ (200.5KB, xls)
   File S2

   The details of differentially expressed genes at +3 d compared with −7
   d.

   (XLS)
   [153]Click here for additional data file.^ (287.5KB, xls)
   File S3

   The details of differentially expressed genes at +3 d compared with −35
   d.

   (XLS)
   [154]Click here for additional data file.^ (399KB, xls)

Acknowledgments