Core circadian clock transcription factor BMAL1 regulates mammary epithelial cell growth, differentiation, and milk component synthesis

The role the mammary epithelial circadian clock plays in gland development and lactation is unknown. We hypothesized that mammary epithelial clocks function to regulate mammogenesis and lactogenesis, and propose the core clock transcription factor BMAL1:CLOCK regulates genes that control mammary epithelial development and milk synthesis. Our objective was to identify transcriptional targets of BMAL1 in undifferentiated (UNDIFF) and lactogen differentiated (DIFF) mammary epithelial cells (HC11) using ChIP-seq. Ensembl gene IDs with the nearest transcriptional start site to ChIP-seq peaks were explored as potential targets, and represented 846 protein coding genes common to UNDIFF and DIFF cells and 2773 unique to DIFF samples. Genes with overlapping peaks between samples (1343) enriched cell-cell adhesion, membrane transporters and lipid metabolism categories. To functionally verify targets, an HC11 line with Bmal1 gene knocked out (BMAL1-KO) using CRISPR-CAS was created. BMAL1-KO cultures had lower cell densities over an eight-day growth curve, which was associated with increased (p<0.05) levels of reactive oxygen species and lower expression of superoxide dismutase 3 (Sod3). RT-qPCR analysis also found lower expression of the putative targets, prolactin receptor (Prlr), Ppara, and beta-casein (Csn2). Findings support our hypothesis and highlight potential importance of clock in mammary development and substrate transport.


Introduction
Circadian clocks set daily rhythms of gene expression to maintain tissue homeostasis and coordinate cellular metabolism [1][2][3]. In mammary tissue, temporal expression patterns of core clock genes change across reproductive states of the female, and appear to be associated with changes in developmental stage of the gland [4,5]. However, the significance of these changes and the role of the molecular clock and core clock genes' functions in mammary development and lactation is currently unknown. Clocks' molecular mechanism generate circadian rhythms through a series of interlocked transcription-translation feedback loops. CLOCK and BMAL1 are at the core of the loop, and as a heterodimer (BMAL1:CLOCK) function as a transcription factor that binds the enhancer box (E-box) regulatory element in promoter regions of genes [6]. Period (PER1, PER2 and PER3) and Cryptochrome (CRY1 and CRY2) genes are transcriptional targets of BMAL1: CLOCK and form the negative arm of the core circadian clock loop, whereby PER and CRY proteins together inhibit CLOCK:BMAL1-mediated transcription [7][8][9][10]. The 24 hr periodicity in activation-repression of elements of the core molecular clock results in corresponding circadian rhythms of gene expression in about 10% of the transcriptome [8].
Although there is some overlap, the rhythmic output genes of clocks differ among tissues, allowing circadian control of function and activity appropriate for each organ. Circadian and cell cycle regulation is coupled to coordinate cell proliferation with tissue function across multiple organs [8,11,12]. This coupling likely also exists in the mammary gland, as temporal analysis of mammary transcriptomes of mature virgin mice [13] and lactating women [14] demonstrated genes that regulate cell growth and differentiation exhibit circadian rhythms of expression. Moreover, the Clock-Δ19 mice, which have a point mutation that results in down regulation of CLOCK-BMAL1 target genes, exhibit decreased lactation competence as evident in lower rates of pup growth and increased rates of neonatal mortality [15][16][17]. The decreased postnatal pup survival was associated with impaired mammary development in Clock-Δ19 dams [18], and isolated mammary epithelial cells from virgin Clock-Δ19 mice had lower stem cell-like properties [13]. Whereas reduction of CLOCK protein in a mammary epithelial cell line with shRNA affected cell growth and decreased expression of factors associated with mammary differentiation and milk synthesis [18].
With this knowledge, we hypothesized that circadian clocks in the mammary gland play an integral role in regulating mammogenesis. In particular, we propose that the transcription factor BMAL1:CLOCK functions to regulate expression of genes that control mammary epithelial development. Moreover, the observation that stoichiometry changed between positive (BMAL:CLOCK) and negative (PER2) elements of core circadian clock components in the mammary gland during the transition from pregnancy to lactation [5] led us to the hypothesize that transcriptional targets of BMAL1:CLOCK change as the gland differentiates and function to regulate milk component synthesis. The objectives of this study were to use ChIP-seq analysis to identify transcriptional targets of BMAL1, and determine if targets of BMAL1 changed upon differentiation of mammary epithelial cells. A line of mammary epithelial cells with the BMAL1 gene knocked out (BMAL1-KO) using CRISPR-CAS9 was created to verify ChIP-seq findings, and study functional effects of loss of gene on growth and differentiation of cells in culture.

Validation of BMAL1 antibody and quality of ChIP-seq data
To test the hypothesis that BMAL1 transcriptional targets changed with state of mammary epithelial differentiation we used ChIP-seq analysis to identify genes potentially regulated by BMAL1 in undifferentiated (UNDIFF) and lactogen hormone differentiated (DIFF) HC11 cultures. Prior to beginning studies, specificity of the ChIP-grade antibody for the BMAL1 protein was confirmed with immunoprecipitation and western blot analysis (S1A Fig). Nanochip analysis of sonicated input DNA indicated that optimal size of sheared fragments was achieved for ChIP-seq analysis (S1B Fig). Evaluation of antibody specificity indicated no difference between mock-ChIP and BMAL1-ChIP samples following RT-qPCR analysis of an exon region of the Magea1_2 sperm specific gene, which is not a BMAL1:CLOCK target. Whereas a 9-fold difference in enrichment was found between RT-qPCR product of BMAL1-ChIP and mock-ChIP for the Per1 promoter region versus the exon region of Magea1_2 (S1C Fig).
The eight samples sequenced for these studies consisted of four pairs of input and ChIP samples with two pairs from UNDIFF and two from DIFF HC11 cultures. Across the four samples the amount of DNA captured by ChIP was 1.03 ± 0.24% of the input DNA. Mapping rate of reads to the mouse genome across the four paired input-ChIP samples averaged 96% (S1 Table). For UNDIFF samples, there was an average of~6,000 peaks with a 2-fold enrichment over input tag count (FDR � 0.01), and over 13,000 peaks were identified across both DIFF samples. Annotation analysis of peaks found similar frequencies of location within the genome across the samples, with approximately 62% of peaks located in intergenic sites and 34% within introns (Fig 1A). Less than 1% of peaks were found in 3' untranslated region, transcriptional termination sites (TSS), and within exons. Functional analysis of protein coding genes nearest transcriptional start sites of ChIPseq peaks. Ensembl protein coding gene IDs with the nearest transcriptional start site to peaks were explored as potential regulatory targets of BMAL1. Analysis of overlapping gene IDs among the samples found 846 common to all four ChIP-seq samples (i.e. UNDIFF and DIFF; Fig 1B; S2 Table). There were 2773 protein coding genes common to both DIFF samples, but not found in UNDIFF samples, and were considered as potential BMAL1 targets distinctly regulated in differentiated versus undifferentiated mammary epithelial cells (S3 Table). To increase confidence in potential as a BMAL1 target, only genes with overlapping peaks between at least two samples were used for downstream analysis. There were 97 common between UNDIFF samples (S4 Table) and 778 common between DIFF samples (S5 Table). Gene targets with overlapping peaks between any two samples (1343; S6 Table), regardless of differentiation state, were used for downstream analysis.
Functional annotation analysis of the 1343 genes with overlapping peaks (Table 1; S7 Table) found approximately one-third enriched the gene ontology membrane. Among the genes in this category were nine glutamate receptors, ten solute carriers, the Wnt receptor Fzd1 and receptors for prolactin, growth hormone, insulin and parathyroid hormone. Within the category GO:0007156~homophilic cell adhesion via plasma membrane adhesion molecules were fourteen cadherin genes and three atypical cadherins. Also, among the BMAL1 targets were genes that encoded multiple ephrin receptors and their ligands, prostaglandin E receptors, mitochondrial fission regulators, cell division regulators, multiple steroid hormone receptors and extracellular matrix proteins and proteases (Table 1).
Ingenuity Pathway Analysis (IPA) tools were used to further explore the data set of potential BMAL1 targets represented by overlapping peaks between samples. The most enriched IPA canonical pathway (S8 Table) was Synaptogenesis Signaling Pathway (p = 1.21E-10; S9 Table). Among the 43 genes that enriched this pathway were fourteen cadherins, ephrins and their receptors, glutamate receptors, the very low density lipoprotein receptor (Vldr), and alpha synuclein (Snca). Another highly enriched canonical pathway was Axonal Guidance Signaling (p-value = 4.67E-6; S10 Table). Genes that enriched this pathway included two round about guidance receptors (Robo1, Robo2) and their ligand (Slit2), semaphorins (Sema3A, Sema5A), Wnt2, Wnt2b, unc-5 netrin receptors (Unc5c, Unc5d) and netrin G (Ntng1). The most significant upstream regulators of BMAL1 (S11 Table) target genes identified by IPA were the glutamate receptor GRIN3A (p-value = 1.05 E-16; S12 Table), the CREB transcription factor (p-value = 2.76 E-14; S13 Table), and the chorionic gonadotropin complex (CG; pvalue = 2.28 E-8; S14 Table) with 32, 72 and 45 molecules in the datasets, respectively. The IPA generated regulator network with the highest consistency score indicated the potential BMAL1 targets positively affected processes involved in formation of cellular protrusions, branching of cells, development of sensory organ, sprouting, size of body, and efflux of lipids, while inhibiting organismal death (S13 Table; S2 Fig) A schematic, which represented approximately 15% of the BMAL1 targets, was generated using IPA my pathway tools by combining genes identified as downstream targets of CREB and estradiol (Fig 2; S14 Table). The schematic illustrates that BMAL1 targets encode proteins across multiple subcellular locations, and the presence of multiple ligand and receptor pairs to include prolactin (Prl8a2) and its receptor (Prlr), Slit2 and Robo, Wnt2 and Fzd1, glial cell derived neurotrophic factor (Gdnf) and its receptor (Gfra1) and insulin (Ins1) and its receptor (Insr).

Analysis of overlap of protein coding genes closest to ChIP-seq peaks with genes exhibiting circadian rhythms of expression in mammary and liver BMAL1 target genes
A secondary approach to verify BMAL1 targets identified with ChIP-seq, was to query for overlap of targets with genes that exhibited circadian rhythms of expression in mammary tissue or identified as BMAL1 targets in hepatic tissue. The set of genes in S2 and S3 Tables were used for this analysis. Genes identified as potential targets in UNDIFF and DIFF samples overlapped with 102 genes that exhibited circadian rhythms of expression in virgin mouse glands [13] and 189 in lactating breast of women [14]. Functional annotation analysis of overlap between mammary circadian rhythm transcriptomes and potential BMAL1 targets categorized thirteen genes in GO:0006974~cellular response to DNA damage stimulus (S17 Table and  Table 2). Nine genes were clustered in GO:0006629~lipid metabolic process including, Fdft1, Gpam, a thromboxane synthase (Tbxas1), and several lipases (Mgll, Plcg2, Plcl2, Pla2g4a).
https://doi.org/10.1371/journal.pone.0248199.g002 Table 2. Representative categories enriched with potential BMAL1 target genes identified using ChIP-seq that overlapped with genes that showed circadian rhythms of expression in mammary glands and BMAL1 targets identified in hepatic tissue.

Term Count % P-value Genes
Overlap with genes that showed circadian rhythm of expression in mammary glands the cell cycle across the eight days ( Fig 4E). FACS screening for cells with less than 2N, an indicator of dead or dying cells, found a higher percent in BMAL1-KO cultures on all eight days (p<0.05; Fig 4F). Despite no difference in the percent of cells in S/G2/M phases, the expression level of Ccnd1, which regulates the transition from G1 to S phase of the cell cycle, was greater (p<0.05) in BMAL1-KO line ( Fig 4G).
Others reported that transgenic mice with BMAL1 gene knocked out (Bmal1 −/− ) had reduced lifespans and displayed symptoms of premature aging associated with increased levels of reactive oxygen species in some tissues [20]. To determine if this was the potential cause of cell loss in BMAL1-KO line, reactive oxygen species (ROS) were measured and found to be significantly higher in BMAL1-KO cultures on day 3 ( Fig 4H). RT-qPCR analysis of the antioxidant Sod3, which was identified as a potential target of BMAL1 in HC11 cells (one UNDIFF and both DIFF cultures, peaks did not overlap), revealed that mRNA levels were significantly lower in BMAL1-KO versus wild-type HC11 cultures ( Fig 4G).
In the DIFF2 sample the serotonin transporter-SERT (Slc6a4) and tryptophan hydroxylase 1 (Tph1), which encodes the protein that catalyzes the first and rate-limiting step in the biosynthesis of serotonin, were identified as potential targets of BMAL1. Lactogen induced differentiation of HC11 cells significantly increased mRNA levels of Tph1 ( Fig 5A) and Sert (Fig 5B), whereas levels were significantly lower in the BMAL1-KO line. RT-qPCR analysis using primers that targeted two sites in the Sert/Slc6a4 promoter region that contained E-box sequences beginning at -42 and -1282 nucleotide bases upstream of the transcriptional start site, found levels 3.5-fold and 2-fold, respectively, higher in UNDIFF BMAL1 ChIP samples than mock-IP samples (Fig 5C). Although only one DIFF sample had a peak in the promoter region of the Sert/Slc6a4 gene (S4F Fig).
Ppara, which regulates lipid metabolism and glucose homeostasis, was also identified as a potential target in both DIFF samples (S4A and S4B Fig) as well as in the liver [19]. RT-qPCR analysis of Ppara expression levels found mRNA depressed in UNDIFF and DIFF BMAL1-KO cultures relative to wild-type HC11 (Fig 6A). Fatty acid synthase (Fasn) was identified as a For this experiment cells were grown to confluence in growth media. Media was changed to lactogen media for 2 hr to synchronize clocks, as we previously have shown this synchronizes clocks in HC11 cells [5]. At completion of 2 hr lactogen treatment (time 0 hr), cells were rinsed with PBS and cultured in growth media for remainder of the experiment. Cells were collected for isolation of total RNA every 4 hr over a 48 hr period beginning at 0 hr. Per2 was measured with RT-qPCR, and levels were expressed relative to mean ΔCT of HC11 across all time points. Cosinor analysis found mesor (-0.12 and -0. Although Csn2, which encodes the milk protein beta-casein, was only identified as a potential target in one of the differentiated samples (DIFF2), levels of mRNA ( Fig 6B) and protein ( Fig 6C) were lower in BMAL1-KO versus wild-type HC11 lines in DIFF cultures. Prolactin regulates Csn2 expression and the prolactin receptor (Prlr; S4C and S4D Fig) was identified as a BMAL1 target. Prlr mRNA expression was reduced in UNDIFF and DIFF BMAL1-KO cultures ( Fig 6D). Lower Prlr levels would affect multiple pathways that stimulate mammary  Values are mean across triplicate samples and two replicate experiments, normalized to express fold-change relative to mean of HC11 ± standard deviation using delta-delta cycle threshold method. Note the y-axis in Ppara is fold change, whereas y-axis for Csn2 is log base 2 of fold change due to large induction. ANOVA and post-hoc Tukey test analysis findings is indicated by differing letter reflecting difference at p<0.05. (c) ELISA quantification of CSN2 protein in UNDIFF and DIFF HC11 (black) and BMAL1-KO (gray) cultures. Values are mean concentration ± standard deviation across triplicate samples. ANOVA and post-hoc Tukey test analysis findings is indicated by differing letter reflecting difference at p<0.05. (d) RT-qPCR analysis of Prlr in UNDIFF and lactogen DIFF HC11 (black) and BMAL1-KO (gray) cultures. Values are mean across triplicate samples, normalized to express fold-change relative to mean of HC11 ± standard deviation using delta-delta cycle threshold method. Note the y-axis is fold change. ANOVA and post-hoc Tukey test analysis findings is indicated by differing letter reflecting difference at p<0.05. (e) Images of two and a half dimensional drip gel cultures of HC11 and BMAL-KO cells taken with phase-contrast microscopy after 7 days of incubation in lactogen media. Cells were plated at 13,000 cells/well.  Values are mean across triplicate samples and two experimental replicates. To calculate relative difference, data were normalized to mean ΔCT of HC11 and fold change difference was determined using the delta-delta CT method; ANOVA and post-hoc Tukey test analysis indicated with differing letter reflecting difference at p<0.05. (c) ChIP-qPCR analysis of Slc6a4 (Sert) promoter region using primers that targeted two sites that contained Ebox sequences beginning at -42 and -1282 nucleotide bases upstream of the transcriptional start site in undifferentiated HC11 cultures; values are mean across four samples, normalized to express fold-change relative to mean CT of mock. A 2-fold difference was considered a positive ChIP.

Discussion
Genes identified in undifferentiated and lactogen differentiated mammary epithelial cells as potential transcriptional targets of BMAL1 support the hypothesis that the circadian clock in the mammary gland regulates mammary development and lactation, and that transcriptional targets of BMAL1 change as mammary epithelial cells differentiate. The greater number of transcriptional targets in differentiated HC11 cultures were likely due to the effects of lactogenic hormones on chromatin access. Higher doses of glucocorticoid treatment are associated with pioneering activities of glucocorticoid receptors. Binding of activated glucocorticoid receptors to chromatin increases the accessibility of regions of the genome previously inaccessible to other transcription factors [22,23]. Prolactin (Prl8a2) and the prolactin receptor (Prlr) were identified as BMAL1 targets. Prolactin regulates Bmal1 expression in mammary epithelial cells [5], and mammary epithelial cells synthesize and secrete prolactin [24]. Thus, in differentiated mammary epithelial cells BMAL1 may play a role in the iterative induction of its own activity through positive loops with prolactin signaling. Upon prolactin binding to its receptor JAK2 is activated. Jak2 was also identified as a BMAL1 target. Activated JAK2 phosphorylates and activates STAT5 transcription factors. Activated STAT5 transcription factors were associated with pioneer-like effects on chromatin in mammary epithelial cells [25]. Thus, the higher number of BMAL1 targets in the lactogen differentiated cultures may have been due to dexamethasone and prolactin treatments increasing the accessibility of BMAL1 to chromosomal regions. Moreover, the positive loop between BMAL1 and prolactin signaling molecules may potentially explain the changes in stoichiometric relationships between positive and negative elements of the mammary epithelial clock transcription-translation feedback loop observed following lactogen induced cellular differentiation [5].
The prediction of CREB as an upstream regulator of BMAL1 targets is consistent with CREB's role as a coactivator of BMAL1:CLOCK transcription factor activity [26,27]. Genes identified as CREB targets encompassed cell adhesion molecules (cadherins, catenin and nectin genes), leptin, lipoprotein lipase, and estrogen receptor gamma. The prediction of the glutamate-regulated ion channel GRIN3A as the most significant upstream regulator of BMAL1 target genes, suggest that, similar to the role that glutamate places in regulation of the master clock in the SCN [28], it potentially entrains circadian rhythms in mammary epithelial cells. Moreover, the identification of chorionic gonadotropin complex as an upstream regulator of BMAL1 target genes may reflect a potential role of the mammary epithelial clock in integration and coordination of the timing of endocrine-paracrine signal-receptivity-response. Among the targets downstream of this complex were receptors for mineralcorticoid (Nr3c2), prolactin (Prlr), prostaglandin E (Ptger2), natriuretic peptide (Npr3), luteinizing hormone/choriogonadotropin (Lhcgr), bone morphogenetic protein (Bmpr1b). In addition, there were receptors for Wnt ligands (Fzd1), ephrins (Epha3 and Epha7), an integrin (Itgb1), and the orphan nuclear receptor Nr5a2, which acts as a metabolic sensor that regulates expression of genes involved in bile acid synthesis, cholesterol homeostasis and triglyceride synthesis [29,30].
Among the BMAL1 targets identified were multiple components of the Wnt signaling pathway (Wnt2, Wnt2b, Fzd1, Rspo2, Ctnna2) and Cdk14, which acts as a cell-cycle regulator of Wnt signaling pathway. The Wnt signaling pathway influences mammary stem cells to elicit fate specification and patterning during embryonic and postnatal development of the gland [31]. Regulation of the Wnt/β-catenin pathway by the circadian clock has been demonstrated in epidermal stem cells [32]. Alterations in expression level of genes in the Wnt pathway may potentially be the cause of compromised mammary epithelial stem cell function observed in Clock-Δ19 mice [13], as well as the premature aging evident in Bmal1-/-mice [33]. The nuclear hormone receptors for glucocorticoids, retinoids, and estrogens also regulate adult stem cell fate [33]. Retinoic acid receptor beta (Rarb), the orphan receptor estrogen-related receptor gamma (Esrrg), which interferes estrogen receptor signaling, and the mineralcorticoid receptor (Nr3c2), which binds glucocorticoids, were found among the potential targets of BMAL1 in mammary epithelial cells. Additional evidence that the clock in mammary epithelial cells may regulate cell fate and developmental patterning was the network generated by IPA software that indicated BMAL1 targets affected branching of cells, development of organs, and body size, while inhibiting cell death. Moreover, numerous BMAL1 targets were involved in stem cell maintenance to include Sox1 and Sox6 [34] and Pard3 and Pard3b genes. Pard3 is necessary for mammary gland morphogenesis, and Pard3b is expressed by multipotent stem cells in the terminal end buds of murine mammary glands, with studies finding ablation of Pard3b resulted in stem cell loss [34].
Also among the BMAL1 targets were a multitude of ion transporters, solute carriers, glutamate transporters and ATP-binding cassette transporters. Synthesis and secretion of milk is dependent on membrane transport systems that move ions and substrates into and out of epithelial cells [35,36]. Movement of ions by ion transporters creates potential differences that enable electrical signaling which regulates cell number, shape, differentiation, and morphogenesis [37]. Glutamine is the primary nitrogen donor for proliferating cells, and alterations in intracellular glutamine abundance was associated with quiescent and proliferative states of mammary epithelial cells [38]. Glutamine also contributes carbons to biosynthetic reactions, and thus mammary clock regulation of glutamate transport may be a means by which the timing of proliferation and metabolic activities of epithelial cells are coordinated. Of the ten ATP binding cassette genes identified as BMAL1 targets seven were the ABCA subfamily type, which function to mediate efflux of cholesterol (Abca1) and other lipids from cells. Others included Abcc12, which mediates toxin efflux, the xenobiotic transporter Abcg2, and Abcb5, which may also function in the efflux of drugs and toxins. Thus, enrichment of pathways and categories related to ion channels, synapse and substrate transports point to the potential role of the mammary clock in affecting the program of mammary development as well as regulating substrate availability and detoxification during lactation.
Functional studies with the BMAL1-KO line found that deletion of Bmal1 did not affect proliferation rate of cells, but rather resulted in greater rates of cell death and lower metabolic activity. Although, Ccnd1 was identified as a potential BMAL1 target, deletion of Bmal1 gene resulted in significantly higher expression of Ccnd1 in BMAL1-KO cultures. Similarly, when abundance of CLOCK protein was decreased in HC11 cells with shRNA, Ccnd1 expression was higher and proliferation rate increased [18]. In our previous study, lower CLOCK abundance appeared to result in the loss of gating of entry and progression through the cell cycle. Elevated Ccnd1 likely reflected this phenomena in culture, and thus loss of gating may also explain elevated Ccnd1 levels in BMAL1-KO cells. Gating the time cells enter the cell cycle by clocks enables the temporal separation of processes incompatible with DNA synthesis [39]. Loss of gating time, coupled with increased levels of reactive oxygen species, related to lower Sod3 expression, in BMAL1-KO cells may have resulted in genotoxic stress, and cell death. There were 31 BMAL1 targets that encoded proteins involved in response to DNA damage, to include Cdkn2aip, Nsmce2, Rad23b, Hus1, Ubr5, Ube2e2, Ube2w. Thus, death of BMAL1-KO cells was also potentially due to decreased ability to repair DNA damaged by reactive oxygen species. Moreover, Ubr5, Ube2e2, Ube2w are involved in the ubiquitination pathway. The ubiquitination pathway targets proteins for degradation and plays a central role in maintaining proteostasis, the cellular balance of protein synthesis and degradation. The aging process is related to changes in proteostatic equilibrium and build-up of misfolded and damaged protein aggregates [40,41], and thus the observation that Bmal1 knockout mice (Bmal1 -/-) exhibited advanced aging [20], may be due to alterations in these genes and the ubiquitination pathway.
Involvement of BMAL1 targets in proteostasis supports that the mammary epithelial clock regulates cellular-tissue homeostasis. During lactation, the serotonergic system plays a central role in regulating mammary epithelial cell homeostasis [42]. Multiple components of the serotonergic system were identified as potential BMAL1 targets in several of the samples. TPH1 and SERT/SLC6A4 are central to serotonergic regulation of mammary epithelial homeostasis through regulation of synthesis and degradation of serotonin, respectively. Our previous analysis found multiple E-box sequences in the promoter region Sert/Slc6A4 and Tph1 genes, and that SERT exhibited circadian rhythms of expression in HC11 cells and lactating sheep mammary [43]. ChIP-qPCR analysis of regions encompassing E-box sequences in the promoter regions of SERT confirmed BMAL1 binding to these sites in HC11 cells. Levels of Sert and Tph1 were also significantly depressed in BMAL1-KO cultures relative to wild-type HC11. The circadian system functions in part to maintain metabolic homeostasis [44], and thus a role for the mammary clock in regulating factors that control epithelial homeostasis during lactation is consistent with this function.
Although not significantly enriched, multiple BMAL1 targets were categorized as lipid metabolic process. Some overlapped with genes that showed circadian rhythms of expression in mammary and those that were confirmed targets of BMAL1 in hepatic tissue. Notable among them were Fdft1, Gpam, Vldlr, Plin5 and Mgll, all of which increase significantly upon secretory activation of the mammary gland [45,46], a process marked by the onset of milk fat synthesis and secretion. FDFT1 (farnesyl-diphosphate farnesyltransferase 1) regulates cholesterol synthesis, GPAM catalyzes the synthesis of glycerolipids, VLDLR and MGLL function in the uptake and hydrolysis of triglycerides, and PLIN5 mediates lipid droplet formation. Together supporting that the mammary epithelial clock plays a central role in the regulation of milk fat synthesis during lactation.
Also central to lipid metabolism is PPARA (peroxisome proliferator-activated receptoralpha). As a nuclear receptor it senses hormonal and nutrient status of the cell and functions to stimulate uptake, utilization, and catabolism of fatty acids [47]. RT-qPCR analysis of Ppara expression showed that levels were significantly depressed in UNDIFF and DIFF cultures of BMAL1-KO cells versus wild-type HC11, and failed to increase following four days of lactogen treatment. These findings are consistent with the knowledge that PPARA is a well characterized target of the BMAL1:CLOCK transcription factor and believed to be central to coordinating cellular energy metabolism with molecular clocks [47]. The interaction of the mammary clock and PPARAmay explain changes observed in milk composition that occurred with night restricted feeding in dairy cattle [48], although exploration and confirmation of a role for PARA regulation of lipids synthesis in the mammary gland is limited [49]. The transport of nutrients also appears to be under the control of the mammary epithelial clock as multiple transporters showed overlap with BMAL1 targets in hepatic tissue to include citrate, amino acid, lactate, pyruvate, and glucose transporters, and thus potentially also affected circadian rhythms of milk component synthesis.
Fatty-acid synthase (FASN) is the enzyme that catalyzes the first committed step in fattyacid biosynthesis, and products of FASN mediated synthesis serve as PPARA ligands [50,51]. Despite FASN being identified as a BMAL1 target in HC11 cells, there was lack of a significant effect of deletion of BMAL1 on Fasn expression. This finding maybe due to the integrated and reciprocal regulation of FASN products and PPARA activation [50][51][52][53].
BMAL1 expression is significantly increased during the transition from pregnancy to lactation in mouse mammary glands, and this increase is modelled in lactogen treatment of HC11 cells (results and [5]). The increase in BMAL1 is likely a direct response to increased prolactin and prolactin induced signaling in lactogen treated cells and at the onset of lactation, as discussed above. The prolactin receptor (Prlr) was identified as a BMAL1 transcriptional target in HC11 cells, and levels of mRNA were significantly depressed in undifferentiated and lactogen differentiated BMAL1-KO cultures. Lower levels of Prlr would affect all pathways regulated by prolactin, which is a key hormone in induction of mammary epithelial cell growth, differentiation and milk component synthesis [54]. The lower expression of the milk protein beta-casein (Csn2) in DIFF BMAL1-KO cultures, and the reduced ability of this line to form acini in culture were likely due, at least in part, to the lower level of prolactin receptors.
Cell adhesion, cell junction and proteolysis were categories highly enriched with BMAL1 targets identified in HC11 cells. Cell-cell and cell-extracellular matrix (ECM) interactions play a central role in mammary morphogenesis and regulation of milk synthesis [43]. BMAL1 targets within cell adhesion included fourteen cadherin molecules and three atypical cadherins. Cadherins make-up the extracellular and transmembrane component of adherens junctions and function as cell-cell adhesion receptors that mediate interactions between adjacent cells. The cytoplasmic component of adherens junctions is a multiprotein complex that links adherens junctions with the actin cytoskeleton [55]. The actin cytoskeleton is hypothesized to function to relay SCN-driven timing cues to gene expression in peripheral tissues [56]. Thus, it is not surprising that the mammary clock would control cellular components that connect to factors that transmit temporal information within cells. The dramatic decreased ability of BMAL1-KO cells to form acini in culture likely reflect loss of activity in cell-cell and cell-ECM adhesion molecules that are critical to the morphogenesis and differentiation of mammary epithelial cells. Cell-cell and cell-ECM junctional complexes sense environmental changes and function to maintain mammary epithelial homeostasis [57]. Changes in mechanical tension affect circadian clocks in mammary glands of virgin mice [13]. Adhesion molecules and proteolytic enzymes can affect the environment of the cells, and it has been hypothesized they mediate circadian timing in the central nervous as they can sense rapid signaling events and transmittransfer signals into broader changes in cell activity with the daily changes [58].

Study limitations
There are several limitations to our studies. When the frequency of peak location of BMAL1 targets found in our study were compared with findings of others, the rate of intergenic sites was higher in our data set and the frequency of peaks found in promoter regions was lower [59]. The frequency of regions of peaks identified were consistent across the four samples analyzed in our study. Therefore, the difference may be tissue or cell type specific, as our studies were conducted using a mammary epithelial cell line, whereas the other group conducted studies using liver tissue from mice. Although ENCODE guidelines indicated that two replicates are adequate, there was variability across the samples [60]. Additionally, to cast a wider net for identification and functional annotation of BMAL1 target genes, the stringency for peak calling was relaxed from the default settings of HOMER, and in doing so the number of false-positives likely increased. To increase the confidence in BMAL1-targets, genes with overlapping peaks were highlighted in the manuscript, the majority of these were the 846 common to all 4 samples When genes identified in our study were compared with datasets of genes that showed circadian rhythms of gene expression in mammary tissue or identified as BMAL1 targets in hepatic tissue, overlapping between datasets was only realized when criteria were further relaxed to common genes identified, but peaks did not necessarily overlap between samples. Similarly, only when this more relaxed selection criterion was used, were BMAL1 (ARNTL) and CLOCK identified as significant upstream regulators predicted by IPA.
Another limitation to the studies described here, is that circadian rhythms of gene expression and BMAL1 binding activity were not measured. Although robust circadian rhythms of core clock genes are evident in UNDIFF HC11 cells, constant exposure to lactogenic hormones, which is required to differentiate cells, results in loss of circadian rhythms of expression of multiple core clock genes [4,5]. Loss of circadian rhythms of core clock genes' expression also occurred in mice in early lactation and resulted in relatively constant levels of BMAL1:CLOCK transcription factor. This dynamic is likely due to the hormonal milieu of the physiological state [5]. Lack of capturing a circadian rhythm of BMAL1 binding may have led to our findings that expected targets such as Per1 were not identified across all four samples (S4G Fig). Additionally, data were completely generated using a mammary cell line, and need to be confirmed and studied in vivo.
Caution must also be used in interpreting comparisons between HC11 and BMAL1-KO lines. An appropriate negative control was not used for studies as off target effects were observed when the scrambled sequence was transfected into HC11 cells. Off target effects are common with CRISPR-CAS [61], and so findings presented may not be specific to loss of BMAL1 function. The assay used to measure reactive oxygen species level also has limitations and may not reflect true levels of reactive oxygen species content [62], and thus caution is warranted in interpreting these data. Finally, for RT-qPCR analysis 18S was used as the single reference gene for calculating relative gene expression using the delta-delta cycle threshold (CT) method [63]. Whereas, MIQE (Minimal Information for Publication of Quantitative Real-Time PCR Experiments) guidelines call for the use of at least two reference genes [64]. Assessment of 18S as a suitable reference gene for the studies conducted found levels did not vary by cell line, state of cellular differentiation nor across timepoint for analysis of circadian rhythms. The primer-probe set also showed good amplification efficiency (3.3 cycle threshold difference between 10-fold dilutions of samples). Although, it is important to note that the kit used for reverse transcription uses poly-T to improve amplification of mRNA over rRNA. As a ribosomal RNA 18S lacks poly-A and so this may have affected the efficiency of the reaction.

Conclusions
Overall, ChIP-seq analysis revealed potential BMAL1 transcriptional targets in mammary epithelial cells. Knowledge of these targets provides insight into the role of circadian clocks in undifferentiated and differentiated states of mammary epithelial cells. BMAL1 transcriptional targets play central roles in pathways that affect stem cell maintenance, cellular detoxification and proteostasis, substrate transport and milk fat synthesis. There was also evidence for the coordination of endocrine-paracrine signals by the mammary epithelial clock, as BMAL1 targets included both the ligand and corresponding receptor as in the case of prolactin, insulin and Wnts with Fzd receptor. Data presented in the manuscript will serve as a good resource to inform future studies aimed at understanding the role of the mammary clock in mammary morphogenesis and lactation as well as provide insight into the link between circadian disruption and breast cancer and poorer milk production. In vivo studies are now needed to understand the roles of mammary epithelial clocks in stage specific development of the gland in vivo, with interaction of stromal tissue and physiological context maintained.
To induce differentiation, cells were plated at 100,000 cells/mL and grown to confluence. At confluence cells were washed 2-times with PBS, and incubated for 48 hr in RPMI media supplemented with 10% serum and Ins (no EGF in media). Undifferentiated cultures (UNDIFF) were harvested for downstream analysis at this stage of culture. For differentiated cultures (DIFF), media was changed to RPMI supplemented with 10% calf serum, dexamethasone (0.1 μM), insulin (5 μg/mL), and prolactin (5 μg/mL; ovine prolactin: L6520-250IU, Sigma-Aldrich), and incubated for 96 hr, with media change every 2-days. After 96 hr in lactogen cocktail, differentiated cultures were collected.
For growth curve analysis, 100,000 cells were plated using growth media in 6-well dishes, and cells from 2-dishes per time point were harvested using trypsin EDTA and counted on days 2, 4 and 8 using a BioRad TC10 Automated cell counter (BioRad Laboratories, Inc). Growth curves were conducted 5-times, and results were presented as mean number of cells each day ± standard deviation. Doubling time was calculated using the tool available here [65].
To measure circadian rhythms of gene expression, cells were grown to 80% confluence in 6 well plates. Media was changed to lactogen media, and cells were incubated for 2 hr at 37˚C to synchronize clocks [5]. After completion of 2 hr lactogen treatment [designated as circadian time 0 (CT0)], three wells of cells of each line were rinsed with PBS and RLT buffer from the QIAGEN RNeasy kit was added to cultures and lysates were collected and stored at -80˚C for subsequent isolation of total RNA. Cells in remaining wells were washed three-times with PBS, and complete growth media was added. Beginning from CT0 and every 4 hr over the next 48 hr, cells were rinsed with PBS, RLT buffer was added, and lysates collected and stored at -80˚C for subsequent isolation of total RNA.

Two-and a half dimensional drip cultures and immunofluorescence
RPMI growth media supplemented with 5% serum, 5 μg/mL prolactin, 5 μg/mL insulin, and 0.375 ng/μL hydrocortisone (cat. no. 65966, BD Biosciences) was used for 2.5-dimensional drip cultures to study acini formation [66]. Briefly, cold Matrigel (Growth Factor Reduced: 356231, Corning Life Sciences; 50 μL of gel/cm3) was used to coat the bottom of each well of 4-well chamber plates (Lab-Tek, Nunc). The 4-well chamber plate was placed in 5% CO 2 at 37˚C for 30 minutes to allow the Matrigel to solidify, and then cells were plated at a density of 13,000 cells/cm 2 . Cells were allowed to attach Matrigel for 15 min and then media with 5% Matrigel was added drop by drop to cover cultures. Media was changed every two days. On day seven of cultures, images of cultures were captured at 200 X magnification (20 X objective) with phase contrast on a Zeiss AxioVert (Oberkochen, Germany) inverted microscope using a Nikon camera.

Creation and screening of CRISPR-CAS BMAL1 knockout HC11 lines
ORIGENE's ARNTL Mouse Gene Knockout Kit (CRISPR; CAT#: KN301604; Rockville, MD US) was used to knockout (remove) the BMAL1 gene from HC11 cells following manufacturer's protocol. The kit contained: KN301604G1, Arntl gRNA vector 1 (gRNA1) in pCas Guide CRISPR vector with the target sequence: GAACCGGAGAGTAGGTCGGT; KN301604G2, Arntl gRNA vector 2 (gRNA2) in pCas Guide CRISPR vector with the target sequence: CAT-GAAGTCGCTGATGGTTG; KN301604-D, donor DNA containing left and right homologous arms and green fluorescent protein-puromycin resistance gene functional cassette; and GE100003 scramble sequence in pCas-Guide vector. HC11 cultures in 6-well dishes were plated overnight and transfected with the gRNA1, gRNA2 or scramble along with the donor cassette using TurboFectin (Origene, cat# TF81001). Forty-eight hours post-transfection, cells were split 1:10. Cells were split 1:10 a total of seven times. Transfected cells were selected for using puromycin at 8 ug/mL. Monoclonal colonies were established by collecting cells using cloning discs soaked in trypsin. Viable colonies were grown to confluence and tested for heterozygous or homozygous knockout of BMAL1 using PCR and western blot analysis.
PCR primers used to confirm genomic integration were as follows for the left integration junction: TACTAATGTAGCCCAGGATGGT (Sense) and TAGGTGCCGAAGTGGTAGAA (Anti-Sense); right integration junction: AATGGAAGGATTGGAGCTACG (Sense) and CTCAATGAT CTGGGATGACTTACA (AntiSense); and donor cassette CAGATGCCGGTGAAGAAAGA (Sense) and GGAATGAGCTGGCCCTTAAT (AntiSense). PCR amplified DNA was Sanger sequenced at Purdue University's Genomics Core.
Preliminary growth curve analysis studies indicated multiple monoclonal lines created using scramble sequences exhibited off target effects on cell growth. Blasting the scramble sequence against the mouse genome found multiple growth regulatory genes near similar sequences. Alterations in coding sequence, etc., due to integration of scramble could have affected these genes. Thus, a negative control was abandoned for all subsequent studies.

Fluorescence activated cell sorting (FACS) analysis of percent of proliferating cells
For FACS analysis of percent proliferating cells across eight-days of culture, cells were plated at a density of 100,000 cells/ml, and duplicate samples of each line were harvested on days 2, 4 and 8. Cells were counted and approximately 1.5 x10 6 cells were pelleted by centrifugation. Cell pellets of were resuspended in 100 μl of PBS, fixed with drop-wise addition of 280 μl of ice-cold 90% ethanol and stored at -20˚C until day of analysis. On day of analysis, fixed cells were pelleted by centrifugation and resuspended in 1ml PBS-0.5% BSA, pelleted again, resuspended in PBS-0.5% BSA with 100 U/ml RNAse (R6513 Sigma-Aldrich), and incubated for 15 min at 37˚C. Five μl of propidium iodide (PI; 500 μg/ml stock solution) was added and cells were incubated 15 minutes at RT. PI labelled cells were analysed by fluorescently activated cell sorting (FACS) using a Beckman Coulter FC500 instrument in the Purdue University's Flow Cytometry and Cell Separation Facility. FACS analysis was repeated 3-times.

MTT cell proliferation and reactive oxygen species (ROS) assays
HC11 and BMAL1-KO cells were plated in 96-well plates at 3,000 cells/well in RPMI 1640-phenol free growth media. For MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assay, on day 2, 4, 6, and 8 after plating (three wells per line per day), media was removed from wells and cells were incubated with MTT from the Cell Proliferation Kit (Sigma Aldrich, catalogue no. 11465007001), and the colorimetric assay was performed following manufacturer's directions. After completion of the assay, images of cells were captured under bright field at 10 X magnification with a Nikon camera attached to a Zeiss inverted microscope, and absorbance was read at 540 nm on the Spark multimode microplate reader (Tecan Trading AG, Switzerland). To assess relative amount of activity per cell on days 2 and 6 of culture, MTT staining per unit area of cells was measured using the ColonyArea Plugin [67] to capture and calculate colony intensity percent. Colony intensity percent is the ratio of the sum of pixel intensities in a region to the sum of all the pixels within the same region of interest multiplied by 255, i.e. assuming highest intensity with full saturation of these pixels. For analysis 100 individual cells were selected across three images and percent intensity of staining was measured. The Reactive Oxygen Species Assay Kit (Abcam catalogue No. ab113851) was used to measure the level of ROS in cell lines on days 3 and 4 (three wells per line per day) after plating cells, with fluorescence immediately measured at Ex/Em 485/535 nm on the Spark multimode microplate reader. Both assays were performed three times.

Protein isolation, enzyme-linked immunosorbent assay (ELISA), western blot and immunoprecipitation (IP) analysis
Protein lysates were isolated from cultures by pouring off media, washing twice with chilled PBS, and harvesting using a scraper and 3 ml of cold PBS. Cells were pelleted by centrifugation, and cell pellets were lysed for 30 min on ice with 600 μl of Cell Extraction Buffer (Invitrogen, supplemented with 1mM of PMSF and 50 μl/ml of Protease Inhibitor Cocktail, Sigma Aldrich), with vortexing at 10 min intervals. Protein lysates were transferred to microcentrifuge tubes and centrifuged at 13,000 rpm for 10 min at 4˚C. Protein concentration was measured using a nanodrop (ThermoFischer) and the Coomassie Plus Protein Assay kit (Pierce Coomassie Plus (Bradford) Protein Assay; Thermo Fisher Scientific) following the manufacturer's protocol. Samples were stored at -80˚C until further analysis. BMAL1 (cat. no. LS-F39618) and beta-casein (CSN2, cat. no. LS-F13103) protein levels in cell lysates were measured with ELISA kits from LSBio following manufacturer's directions. Since number of cells and protein concentration changed with state of differentiation of cells, data were expressed as micrograms of CSN2 or BMAL1 protein per mg of protein.
For western blot analysis, 100 μg of protein were loaded per lane and electrophoresed on a 10% TGX precast SDS PAGE gel from Bio-Rad. Protein was transferred onto a nitrocellulose membrane, and membranes were blocked using either 5% BSA (Bovine Serum Albumin), and probed for BMAL1, CSN2 and beta-actin proteins using anti-BMAL1 (ab3350, dilution 1:1,000), anti-CSN2 (LS-C373659, Life Spam, dilution 1:200) and Anti-beta-actin (AbCam ab8227, dilution 1:10,000) antibodies, respectively. Blots were washed and then incubated with secondary antibody (ab97051; 1:5000). Membranes were washed and incubated with the detection reagent Clarity Western ECL Substrate (Bio-Rad). Blots were imaged using the ChemiDoc MP system (BioRad). ImageJ was used for densitometric analysis. Density of test protein (BMAL1 and CSN2) were divided by density of beta-actin (BA) band. To enable statistical analysis across multiple gels run to measure CSN2 protein levels, ratio was normalized to relative expression of HC11 UNDIFF.
For IP analysis of specificity of ChIP grade antibody BMAL1 (Abcam; ab3350), rabbit polyclonal antibody to BMAL1 (Abcam; ab3350) was added to 200 μg protein in 200 μl cold PBST (PBS pH 7.4 with 0.02% Tween-20), and rotated overnight at 4˚C. The antibody protein mixture was added to Dynabeads protein A for immunoprecipitation (ThermoFischer) and pipetted gently to resuspend. Samples were rotated for 2 hr at 4˚C. Tubes were placed on magnetic rack to pellet beads, and washed three times with 200 μl cold PBS. Then 15 μl PBS and 15 μl 2X Laemmli sample buffer with 2-mercaptoethanol added per manufacturer's instructions, and gently pipetted to resuspend. Samples were heat at 100˚C for 5 min and placed on magnetic rack to pellet beads. Supernatant was removed and separated using SDS-Page gel for analysis by western blot. Mouse monoclonal BMAL1 antibody (Santa Cruz Biotechnology, Inc., Dallas, TX US; cat. no. sc-365645) was used for western blot analysis.

RNA isolation and real-time quantitative PCR analysis (RT-qPCR)
RNA was collected from cells and isolated using Qiagen's RNeasy kit with an on-column DNase treatment, and quantity was measured using Nanodrop. RIN scores of total RNA following nanochip analysis on an Agilent 2100 Bioanalyzer were found to vary from 8.0-9.0 across all samples. Promega's GoScript Reverse Transcriptase kit was used to reverse transcribe 500 ng of total RNA into cDNA following manufacturer's protocol. Gene expression was analysed with TaqMan Assays on Demand using 2✕ TaqMan Gene Expression Master Mix (Life Technologies) using a 1:10 dilution of the cDNA product. The CFX Connect Real-Time PCR Detection System (BioRad) was used to run RT-qPCR, which was initiated with 2 min incubation at 95˚C and then 40 cycles of 95˚C for 15 sec and 60˚C for 1min, and standby was set at 4˚C. Multiple genes (beta-actin, beta microglobulin and 18S) were screened as reference genes (housekeeping gene) for calculating relative expression using the 2 -ΔΔ CT method [63]; 18S was chosen as the reference gene based on its levels staying steady across time, differentiation state and genotype. Mouse-specific Assay on Demand TaqMan assays (Life Technologies) used were Per2 (Mm00478099_m1), Sod3 (Mm01213380_s1), Csn2 (Mm04207885_m1), Tph1 (Mm01202614_m1), Slc6a4 (SERT, Mm00439391_m1), Ccnd1 (Mm00432359_m1), Ppara (Mm00440939_m1), Prlr (Mm04336676_m1), and 18S. To calculate relative expression, mean ΔCT (target gene-reference gene) of wild-type HC11 UNDIFF cultures was used as normalizer for relative expression. For temporal analysis of Per2 expression across 48 hr sampling mean ΔCT of HC11 cultures across all time was used as the normalizer for relative expression. Data are presented as mean fold-change ± standard deviation of relative expression levels, or mean of log base two-fold change ± standard deviation for Csn2 expression.

ChIP assay and RTq-PCR verification of known targets
Use of two independent biological replicates for ChIP-seq analysis to identify transcriptional targets is consistent with ENCODE (Encyclopaedia of DNA Elements) guidelines [60]. For ChIP-seq studies, we performed four independent ChIP experiments, two with undifferentiated cultures and two for differentiated cultures, to enable analysis of the effect of differentiation on BMAL1 targets. Briefly, undifferentiated and differentiated HC11 cultures in 100 mm plates were treated with 1% formaldehyde to cross-link proteins to chromatin by adding 625 μl fresh 16% formaldehyde directly to the dishes containing cells and 10 mL media, and incubated at RT on a shaking platform for 7 min. To quench the formaldehyde, 500 μl of 2.5M glycine was added, and incubated at RT on the shaking platform for 5 min. Cells were washed twice with ice cold PBS, and removed from culture dish using a cell scraper. Cells were pelleted by centrifugation and stored at the -80˚C. Samples were collected for ChIP assays in four replicate experiments (sample IDs UNDIFF1, UNDIFF2, DIFF1, DIFF2).
Nuclei were collected by resuspending each fixed-cell pellet in 10 mL of Rinse 1 (50mM Hepes pH 8.0, 140mM NaCl, 1mM EDTA, 10% glycerol, .5% NP40, .25% Triton x100), and incubating 10 min on ice. Samples were pelleted by centrifugation 1,200 X g at 4˚C for 5 min, and then resuspend in 10 mL CiA NP-Rinse 2 (10mM Tris pH 8.0, 1mM EDTA, 0.5mM EGTA, 200mM NaCl), and then centrifuged at 1,200 X g and 4˚C for 5 min. Pellets were washed with 5 ml of CiA Covaris Shearing Buffer (0.1% SDS, 1mM EDTA pH 8.0, 10mM Tris HC1 pH 8.0) twice with centrifugation at 1,200 X g and 4˚C for 3 min after each wash. For sonication, pellets from approximately 3 million cells each were resuspended in 130 μl of CiA Covaris Shearing Buffer with 10% protease inhibitor cocktail. Pellets were sonicated in a Covaris E210 instrument (Covaris, Inc., Woburn, MA) using the following parameters: Duty Factor 5%, Peak power 105W, Cycles per burst 200, 10 min. Following sonication sheared chromatin lysate was transferred to a new centrifuged at 20,000 X g and 4˚C for 15 min. Supernatant (sheared chromatin stock) was divided into three aliquots for input DNA, mock-immunoprecipitation (IP) and IP and stored at -80˚C (-20C), after removal of an aliquot for analysis of product size in a 2% agarose gel and on the Agilent Bioanalyzer Nanochip.
For immunoprecipitation step, cross-linked chromatin using anti-BMAL1 antibody and mock-IP aliquots were thawed on ice. Samples were precleared by adding 50 μL resuspended A beads (DynaBeads Protein A, Invitrogen cat. # 10001D), and incubating at 4˚C for 1.5 hr with rotation. Tubes were placed on a magnetic rack and liquid was transferred to a new microfuge tube. One-quarter of the volume of 5X IP Buffer (250mM Hepes/KOH pH 7.5, 1.5M NaCl, 5mM EDTA, 5% TritonX 100, 0.5% DOC [4-chloro-2,5-dimethoxy-amphetamine], 0.5% SDS) was added, and 10 μL of BMAL1 antibody (Abcam ab 3350) was added to IP tube, and nothing was added to mock-IP tube. Samples were rotated overnight at 4˚C, and then 50 μL of resuspended A beads were added to each tube, and rotated for 3 hr at 4˚C. Magnetic rack was used to remove A beads. Supernatant was collected transferred to a new tube, and washed with 1mL IP Buffer twice with rotation at RT for 3 min, and collection on magnetic rack following wash. Supernatants were then washed with 1 mL of DOC buffer (10 mM Tris, pH8.0; 0.25 M l (iCl; 0.5% NP40, 0.5% DOC, 1 mM EDTA), followed by 1 mL of Tris EDTA (TE) buffer pH 8. Supernatant was removed and 150 μL elution buffer was added and samples were rotated at RT for 20 min. Supernatant was transferred to a new tube and store at -80˚C until DNA isolation.
Chromatin was isolated from supernatant by adding 3 volumes of TE with 1% SDS and RNAse A (Sigma R6513), and vortexing to mix. Samples were incubated at 37˚C for 30 min, and proteinase K was added, vortexed again to mix, and cross-linking was reversed by incubating at 55˚C for~2.5 hr, and then overnight at 65˚C in PCR tubes in the thermocycler. DNA was extracted the next day by adding 230 μL TE and 330 μL phenol/chloroform, vortexing, followed by centrifugation for 1 min. Aqueous layer was transferred into a new microfuge tube and 10 μL of commercial glycogen (ThermoScientific R0551) was added with 30 μL 3M NaCl. One volume of 100% ethanol was added and samples vortexed to mix. Samples were incubated for two days at -20˚C, and then centrifuged 1 hr at 4˚C to pellet DNA. Pellet was washed twice by removing the supernatant and adding 500 μl 70% EtOH. Pellet was resuspended in 25 μL TE and concentration was measured with Qubit fluorometer (Thermo Fisher Scientific) and a Bioanalyzer 2100 (Agilent), and samples stored at -80˚C until sequencing.
To determine specificity BMAL1 antibody, RT-qPCR analysis of ChIP product, input and mock-IP samples were performed using primers designed to target the promoter region of Per1, a transcriptional target of BMAL1: CLOCK, and the exon region of a mouse sperm gene Magea1_2, which was not expected to be a target of BMAL1 binding. The promoter region of Per1 containing the proximal E-box enhancer was amplified with the following primer set: forward, 5-CCTCCCTGAAAAGGGGTA-3; reverse, 5-GGATCTCTTCCTGGCATCTG-3. Primers used to amplify MAGEA1_2, with forward: 5-GCCTCTGAGTGCTTGAAGAT-3, and reverse: 5-CAGGGCAGTGACAAGGATATAG-3. Primers were designed to promoter regions of SLC6A4 (aka SERT) to encompass an E-box most proximal (beginning at -42 nucleotides) to the transition start site forward: 5-CTCCAGCTGCGGTAGCAGA-3 and reverse: 5-ATTTGTACTTG CGGCCC-3, and more distal (beginning at -1282 nucleotides) Forward: 5-GGAGTTACAGGC ACGGAAG-3 and reverse 5-GCCTGGCCATTCCATGA-3. Triplicate samples were measured using SYBR green real-time PCR master mix (ThermoFischer Scientific). The analysis was conducted using the CFX Connect Realtime PCR Detection System (Biorad) with the following conditions for reactions were 1 cycle at 50˚C for 2 min; 1 cycle at 95˚C for 2 min, 55 cycles at 95˚C for 15 sec, gradient 58-60˚C for 30 sec and 72˚C for 1min; and 1 cycle at 95˚C for 10 sec, melt curve 65˚C-95˚C, increment 0.5˚C for 5 sec. PCR reaction efficiencies were 98% at annealing temperature 58˚C for Per1 and Sert, and at 60˚C the efficiency for Magea1 was 96%.
To calculate relative amount of input target brought down by ChIP with BMAL1 antibody the following equation was use: input sample-C T of ChIP sample (https://www.thermofisher. com/us/en/home/life-science/epigenetics-noncoding-rna-research/chromatin-remodeling/ chromatin-immunoprecipitation-chip/chip-analysis.html). A positive ChIP was defined as at least 2-fold greater than mock-IP sample [68].

Sequencing of DNA, peak identification, and functional annotation analysis
Input and ChIP sample pairs were sequenced together by paired ends reads on the Illumina Novaseq 6000 (San Diego, CA) following the manufacturer's protocols. Where Input refers to the DNA isolated from samples without ChIP assay. Prior to sequencing libraries were prepared using NEXTFLEX Rapid DNA-Seq Library Prep Kit for Illumina Platforms (PerkinElmer, Waltham, MA) according to manufacturer's protocol. The number of reads for three of the input ChIP sample pairs (both UNDIFF samples and one DIFF sample) averaged approximately 100 million reads. The remaining DIFF input-ChIP pair had over 1 billion reads. To account for discrepancy in depth of sequencing, the reads were divided into five groups, and one group of input-ChIP sample pair was used for analysis. Sequence quality was assessed using FastQC (v 0.11.7) and quality trimming was done using Fastx toolkit to remove bases with Phred33 score of less than 30. The reads passing sequence quality criteria and with a length of at least 50 bases were retained for downstream analysis. The quality trimmed reads were mapped against the reference genome (Mus_musculus.GRCm38.chr.fa) using bowtie2 (Version 2.2.9) [69]. ChIP-seq data were submitted to Gene Expression Omnibus (GEO accession number GSE154937).
Peaks were called from mapped files (bam files) using 'callpeaks' script of HOMER (Version 4.10) [70] to detect transcription factor associated peaks using default settings except for the following parameters: fold enrichment over input tag count (F) was set at 2, fold enrichment limit of expected unique tag positions (C) was set at 2 and fold enrichment over local tag count (L) was set at 2 and FDR was set at <0.01. Quality control analyses were run from files generated in read tag directories to include tag count distribution, autocorrelation analysis, genomic nucleotide frequency relative to read positions, and fragment GC % distribution, met criteria for further analysis. Picard Mark Duplicates was performed on each sample (S20 Table). Common nearest gene associated with peaks between UNDIFF1 and UNDIFF2, DIFF1 and DIFF2 were generated using Venny 2.1 [71].
Overlap of BMAL1 targets in HC11 cells with genes that exhibit circadian rhythms of gene expression in virgin mouse mammary glands [13] and in human breasts during lactation [14] was investigated by obtaining supplemental data from these manuscripts. Similarly, overlap of HC11 BMAL1 targets with BMAL1 transcriptional targets identified in mouse liver was also investigated by obtaining supplemental data this manuscript [19]. For comparative analysis, gene symbols were converted to mouse ENSEMBL IDs using the gene conversion tool available in DAVID, and then common genes across data sets were identified using Venny 2.1 [71].
Functional annotation analysis was performed using Database for Annotation, Visualization, and Integrated Discovery-DAVID Bioinformatic Resources 6.8 [72,73], and Ingenuity Pathway Analysis (IPA; QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ ingenuitypathway-analysis). To aid in visualization of IPA generated Figs, genes were assigned a value of 5 and were indicated as red in Fig of Regulator Networks. The network with highest consistency score was used, with the higher scores reflecting paths between target gene and output function consistent with the predicted state of the regulator-here all indicated as upregulated-based on the literature. IPA predicted upstream regulators were removed from the generated network, as BMAL1 was the upstream regulator (S2 Fig). The My Pathway tools in IPA was used for visualization of a subset of BMAL1 target genes. GeneCards [74] was used to query information on function of genes of interest.

Statistical analysis
Statistical analysis was performed using SPSS software (IBM SPSS, v.26). General linear model was used to analyze whether cell line (HC11 or BMAL1-KO) or state of differentiation (DIFF or UNDIFF) affected variables. For time series analysis such as growth curve, MTT assay, ROS assay, and FACS analysis repeated measures was added to model and effect of day and line were analyzed. A Tukey's post-hoc test was used for pairwise comparisons. Significance was considered at P � 0.05. Cosine fit analysis of 24 hr rhythms of gene expression was performed with the cosinor package in R (RStudio 1.1.453, Boston, MA). Mesor, amplitude, acrophase R 2 , and p-value were outputs of the package algorithm.
Supporting information S1 Fig. (a) Western blot analysis of immunoprecipitation (IP) of BMAL1 protein in HC11 cell protein lysate using rabbit polyclonal antibody to BMAL1 (ab3350, ChIP grade, 2 μg per IP). Lane 1 = Precision Plus Protein Standard Ladder; Lane 2 = 100 μg, Lane 3 = 150 μg, Lane 4 = 200 μg, Lane 5 = 200 μg of lysate precleared and incubated with beads but no Ab (-), Lane 6 = an aliquot of supernatant pulled off of lane 5 sample before wash steps that precede elution. western blot was performed using the mouse monoclonal antibody to BMAL 1 (sc-373955 @ 1:750 primary antibody concentration) for visualization. (b) Electropherogram analysis of input DNA used for ChIP-seq shows ideal size for next generation sequencing (seq). (c) Evaluation of antibody specificity indicated no difference between mock-ChIP and BMAL1-ChIP samples in the cycle threshold values following RT-qPCR analysis of an exon region of the Magea1_2 sperm specific gene, which is not a BMAL1:CLOCK target. Whereas a 9-fold difference in enrichment was found between RT-qPCR analysis of BMAL1-ChIP and mock-ChIP for the Per1 promoter region versus the exon region of Magea1_2. A positive ChIP was defined as at least 2-fold greater than mock-IP sample. Different letters indicate significant difference at p<0.05. (TIF) S2 Fig. Ingenuity Pathway Analysis generated network with greatest consistency score reflecting relationships between BMAL1 targets and predicted downstream effects. Gene names and symbols are defined in S15 Table. The predicted upstream regulators were removed, as analysis assumes this is BMAL1 as the transcriptional regulator of genes in red and combined the downstream effects of inhibiting organismal death, and stimulating efflux of lipids, cellular protrusions, branching of cells, development of sensory organ, sprouting and organ size. Values are mean technical replicates within experiment, normalized to express foldchange relative to mean of HC11 ± standard deviation using delta-delta cycle threshold method. Temporal analysis of Fasn expression in WT HC11 (solid line) and BMAL1-KO (dashed line) cultures. For this experiment cells were grown to confluence in growth media. Media was changed to lactogen media for 2 hr to synchronize clocks. At completion of 2 hr lactogen treatment (time 0 hr), cells were rinsed with PBS and cultured in growth media for remainder of the experiment. Cells were collected for isolation of total RNA every 4 hr over a 48 hr period beginning at 0 hr. Fasn was measured with RT-qPCR, and levels were expressed relative to mean ΔCTof HC11 across all time points. Cosinor analysis found mesor (0.01 and 0.97), amplitude (0.21 and 0.12), acrophase (-9.14 and -4.44), R 2 (0.26 and 0.19) and p-value (0.22 and 0.33) of fit to a 24 hr rhythm, respectively, for HC11 and BMAL1-KO lines. (TIF) S1 Table. Total reads, quality reads and overall mapping rate of ChIP and Input samples.  Table. Ingenuity Pathway Analysis generated network with greatest consistency score reflecting relationships between BMAL1 targets and predicted downstream effects. (XLSX) S16 Table. Symbols and corresponding genes names in schematic (Fig 2) generated with Ingenuity Pathway Analysis using my pathway tools. (XLSX) S17 Table. Overlap of BMAL1-Chip seq gene targets that overlap with genes that showed circadian rhythms of gene expression in virgin mouse glands and in lactating breast of women: Gene ID lists and functional annotation analysis. (XLSX) S18 Table. Overlap of BMAL1-Chip seq gene targets with genes that were identified as targets in liver tissue: Gene ID lists and functional annotation analysis. (XLSX) S19