Transcriptional Regulation of Autophagy Genes via Stage-Specific Activation of CEBPB and PPARG during Adipogenesis: A Systematic Study Using Public Gene Expression and Transcription Factor Binding Datasets

Autophagy is the cell self-eating mechanism to maintain cell homeostasis by removing damaged intracellular proteins or organelles. It has also been implicated in the development and differentiation of various cell types including the adipocyte. Several links between adipogenic transcription factors and key autophagy genes has been suggested. In this study, we tried to model the gene expression and their transcriptional regulation during the adipocyte differentiation using high-throughput sequencing datasets of the 3T3-L1 cell model. We applied the gene expression and co-expression analysis to all and the subset of autophagy genes to study the binding, and occupancy patterns of adipogenic factors, co-factors and histone modifications on key autophagy genes. We also analyzed the gene expression of key autophagy genes under different transcription factor knockdown adipocyte cells. We found that a significant percent of the variance in the autophagy gene expression is explained by the differentiation stage of the cell. Adipogenic master regulators, such as CEBPB and PPARG target key autophagy genes directly. In addition, the same factor may also control autophagy gene expression indirectly through autophagy transcription factors such as FOXO1, TFEB or XBP1. Finally, the binding of adipogenic factors is associated with certain patterns of co-factors binding that might modulate the functions. Some of the findings were further confirmed under the knockdown of the adipogenic factors in the differentiating adipocytes. In conclusion, autophagy genes are regulated as part of the transcriptional programs through adipogenic factors either directly or indirectly through autophagy transcription factors during adipogenesis.


Introduction
Autophagy, a self-eating process, is involved in the development and the differentiation of various cell types including the adipocyte. Autophagy contributes to the adipocyte differentiation by reorganizing the intracellular components, forming the fat droplets and/or providing the energy source for the development of the phenotype [1]. Therefore, the knockdown of essential autophagy genes prevents the adipocyte differentiation and intervenes in the lipid accumulation within adipogenic cells [2]. The adipocyte differentiation is a well-understood process. A well defined transcriptional program is initiated and maintained by specific transcription factors. The 3T3-L1 mouse fibroblast is a key model for studying this process. When the pre-adipocyte is induced with a chemical cocktail containing MDI (1-Methyl-3-isobutylxanthine, Dexamethasone, and Insulin), the cell witnesses morphological and metabolic changes and is finally transformed into mature adipocyte [3].

RNA-Seq Dataset
ChIP Several studies reported functional links between the adipogenic factors and key autophagy genes. For example, the transcription factor CCAAT enhancer-binding protein beta (CEBPB) controls the expression of the autophagy-related 4b cysteine peptidase (Atg4b) [4]. A recent study suggested that peroxisome proliferator-activated receptor gamma (PPARG) regulates autophagy through NEDD4 E3 ubiquitin protein ligase (NEDD4) gene expression in thiazolidinediones-treated HepG2 cells [5]. Forkhead box o1 (FOXO1) regulates autophagy induction and it also negatively regulates Pparg, the master adipogenic factor [6,7]. Indeed, the inhibition of FOXO1 is thus necessary for the differentiation to proceed [7,8]. In a previous study from our laboratory, we reported the progressive activation of the Foxo1 gene during the course of differentiation, to reach the highest expression in the mature adipocytes [9]. In addition, we identified a large number of differentially expressed autophagy genes between differentiated and undifferentiated adipocytes. The changes in gene expression were clustered into several groups that corresponded to the known stages of adipocyte differentiation. Therefore, the regulation of autophagy during the course of the differentiation seems to be more dynamic and stage-dependant. Only a few examples of direct and indirect regulation of autophagy genes have been studied.
We hypothesize that the autophagy process may be regulated by the same adipogenic factors either directly or indirectly through specific autophagy transcription factors. To investigate this hypothesis, we modeled the gene expression and transcription factors and co-factors binding at different stages of adipocyte maturation ( Figure 1). Then, we studied the direct regulatory links between the adipogenic transcription factors and key autophagy genes as well as specific autophagy transcription factors. We were able to quantify the changes in the global expression of autophagy genes as the pre-adipocytes progressed toward the maturation stage. We considered changes at both the individual gene and the gene set levels. Then, we investigated the links between the master adipogenic transcription factors and the observed changes in gene expression. We tested whether the established links hold in the knockdown condition of the adipogenic factors. In addition, the role of transcription co-factors and histone modifications in modulating the function of the transcription factor was discussed.

Autophagy Genes are Regulated as Part of the Transcriptional Program of the Adipocyte Differentiation
First, to model the changes in gene expression during the course of the adipocyte differentiation, we used datasets of RNA-Seq 3T3-L1 cells induced with the MDI cocktail and sampled at different time points (−48 to 240 h) ( Figure 1). In agreement with the previous literature, the differentiation course was divided into non, early and late differentiation stages. This classification could explain the variance in expression (>50%) of all genes in multidimensional scaling (MDS) analysis ( Figure 2 and Table S1). Moreover, the variance in the expression of subsets of genes of interest was also explained by the stage of differentiation ( Figure 2). Significant amounts of variance (85% and 87%) in the expression of autophagy genes and autophagy transcription factor genes were explained by the first two dimensions of the MDS (Table S1). Together, the three-stage classification of differentiating cells explains significant amounts of the variance in gene expression. The subset of autophagy genes exhibits temporal changes that correspond to the differentiation course.
To characterize the changes in gene expression in the adipocytes in response to MDI induction, we applied the differential gene expression analysis to all genes after removing low-quality samples and low expressed genes. Significant changes in the global expression were observed when comparing early-differentiated samples to the non-differentiated and these changes were more pronounced in the late-differentiated samples ( Figure S1a). The changes are illustrated by the induction of adipogenesis (>1.5 log 2 fold-change; FDR < 0.2) and lipogenesis (>0.5 log 2 fold-change; FDR < 0.2) markers such as Pparg, CCAAT enhancer-binding protein alpha/beta (Cebpa/b), ATP citrate lyase (Acly), fatty acid synthase (Fasn) and lipoprotein lipase (Lpl) ( Figure S3 and Table S2). Moreover, several lipid metabolism-related gene sets such as lipid droplet and lipid storage were enriched in the early-and late-differentiating adipocytes compared to the pre-adipocytes (Table S3).
The same pattern was observed in the expression of the subset of autophagy genes ( Figure S1b). Key autophagy gene markers such as beclin 1 (Becn1) was up-regulated (0.31 and 0.44 log 2 fold-change; FDR < 0.2) in both early-and late-differentiated cells. Other marker's expression such as microtubule associated protein 1 light chain 3 beta (Map1lc3b) and sequestosome 1 (Sqstm1) was recovered (>0.5 log 2 fold-change; FDR < 0.2) in the late stages compared to pre-adipocytes ( Figure S3 and Table S2).
In addition to changes in the expression of individual gene markers, qualitative changes in subsets of the autophagy genes were observed. In particular, the subset of regulators of macroautophagy and autophagosome maturation were enriched between stages of differentiation (Table S3).

Adipogenic Transcription Regulators Correlate with the Expression of Autophagy Genes
The genes encoding for the adipogenic transcription factors CEBPB and PPARG were one to two log 2 folds up-regulated in the early stage of differentiation (Table 1). In the case of Pparg, the upward expression trend continued to the late stage. Known co-factors such as the mediator complex subunit 1 (Med1) and retinoid X receptor gamma (Rxrg) were also up-regulated in differentiating compared to the non-differentiated cells (Table 1). By dividing autophagy genes into up-and down-regulated groups, we were able to show the potential direction of the regulation if indeed some of these co-expression relations are functional. CEBPB was co-expressed (PCC > 0.25) with down-regulated autophagy genes only in the early stage of differentiation ( Figure 3). By contrast, Pparg was co-expressed (PCC > 0.5) with the regulated genes in both directions and in both early and late stages ( Figure 3). Moreover, the co-factors Rxrg and Med1 showed a co-expression pattern with autophagy genes that is identical or opposite to that of Pparg in at least one of the differentiation stages ( Figure 3).

Master Adipocyte Regulators CEBPB and PPARG Directly Target Key Autophagy Genes
Key autophagy genes such as Map1lc3b and Sqstm1 were briefly down-regulated in non-induced and early differentiated adipocytes. The expression pattern was gradually reversed (0.45 and 0.95 log 2 folds) over time and toward the late differentiation stage. Becn1 was up-regulated (0.31 and 0.44 log 2 fold-change; FDR < 0.2) in both early-and late-differentiated cells ( Figure 4a). The expression of the two adipogenic factors negatively correlates (PCC −0.2 to −0.5) with that of the three autophagy genes. Therefore, they may act as negative regulators or they would be de-recruited away from the regulatory regions of these genes at some point in the differentiation course. Only in the case of Becn1, the correlation is reversed with Cebpb and Pparg in the early and late-differentiation stages ( Figure S4a). Based on this information only, the two factors can work on Becn1 cooperatively or consecutively at their most active stages.
To quantify the active binding of transcription factors on autophagy genes and its change over time, we performed differential peak binding to compare the binding in differentiating sample compared to the non-induced samples ( Figure S2). CEBPB displays a similar active binding on Becn1 and Map1lc3b (Figure 4b). At least one peak around the promoter region of the two genes is observed in non-differentiated cells (Figure 4d). These peaks were less enriched (<−1 log 2 fold-change) in the early differentiation stage only to regain the enrichment to a similar expression in the late stage (Table S4). PPARG show an opposing binding pattern on Becn1 and Sqstm1 (Figure 4b). The change in enrichment is strongest in Sqstm1 where the binding peak in early differentiated cells decreases (−2.75 log 2 fold-change) in the late stage ( Figure 4c and Table S4). The observed patterns suggest dynamic binding which might be a result of the availability or the distribution of the proteins over time. Unlike CEBPB, the binding of PPARG on the key autophagy genes may also occur in fibroblasts and macrophages ( Figure S5a,b).

CEBPB and PPARG Target Autophagy Transcription Factors Genes
With the exception of transformation-related protein 53 (Trp53), several autophagy transcription factor genes were found to be induced in one or both differentiation stages (Table S5). Foxo1, X-box binding protein 1 (Xbp1) and transcription factor EB (Tfeb) were up-regulated (>1.5 log 2 fold-change) overtime to reach the highest expression toward the end of the differentiation process ( Figure 5a). Both adipogenic factors showed active binding patterns on at least one of the three autophagy transcription factor genes ( Figure 4b). The strongest binding activity around the promoter areas of the targets was that of CEBPB on Xbp1 and that of PPARG on Tfeb (Figure 4c,d and Table S4). Both cases were accompanied by positive co-expression of the corresponding genes ( Figure S4b). The two factors seem to also have some binding activity around the promoter of Foxo1 (Figure 4c,d). Most of these binding activities are probably tissue-specific since little to none was observed in other tissue types ( Figure S5).

Auto-Regulation and Transcriptional Feedback Loops May Affect the Regulation of Autophagy
Although not directly related to the regulation of autophagy genes, the adipogenic transcription factors may target their own coding genes to form auto-regulation circuits or a feedback loops (Figure 6b). CEBPB targets the promoter region of its own gene (>1 log 2 fold-change) and that of Pparg (<−0.7 log 2 fold-change) (Figure 6d). The enrichment of the binding peak of CEBPB on its own promoter in the late stage of the differentiation is accompanied by a decrease in the gene expression suggesting a negative auto-regulation (Figure 6a). Together with the binding of PPARG on the promoter of Cebpb (0.56 log 2 fold-change) they form a feedback loop (Figure 6c). The binding of PPARG on Cebpb and the reverse is similarly observed in macrophages ( Figure S5).

Co-Factors and Histone Modifications Modulate the Functions of Adipogenic Transcription Factors
In addition to the changes in binding patterns and intensity of the adipogenic transcription factors, the occupancy of the adipogenic factors is associated with that of the co-factors and the histone modification markers. Specifically, PPARG associates with MED1 and RXRG in all stages of differentiation and all genomic binding regions ( Figure S7a,b). Similarly, PPARG correlates with H3k27ac and H3k4me1 in late differentiation stages and at different genomic locations ( Figure S7c,d). This is not exactly the case for CEBPB where the correlations with the co-factors and histone markers are less at the late-differentiation stage. The same pattern is confirmed by the correlation in occupancy between the transcription factors and other co-factors and histone modifications in autophagy genes. PPARG associates the most with MED1 and RXRG especially in late differentiated cells and at the promoter regions of this subset of genes (Figures 7 and S7). On the other hand, CEBPB is closest to the co-factor E1A binding protein p300 (EP300) at the same stage and location (Figures 7 and S7).
Functionally, MED1 seems to be the most active co-factor of the three. The change in the occupancy of the co-factor mimics that of PPARG on the latter's autophagy targets in late-differentiating cells ( Figure 8a) and CEBPB targets in early-differentiating cells (Figure 8b). In contrast, EP300 and RXRG occupancy changes only correlate with that of CEBPB or PPARG respectively (Figure 8a,b). The histone modification markers that had the strongest correlation with respect to the change in tags on autophagy genes were H3K27ac on CEBPB targets in the early differentiation stage and H3K4me1 on PPARG targets in the late-differentiation stage (Figure 8c,d). . Correlation in occupancy between adipogenic transcription factors and histone markers on autophagy genes at different stages. The total number of reads in peaks of each ChIP-Seq sample (n = 84) was used to represent the factor/histone marker occupancy. Pearson's correlation coefficient (PCC) was calculated for each pair of samples. The PCC value is between 0 (white) and 1 (blue). Samples were categorized by differentiation stage; non, early or late in that order.

Perturbing the Adipogenic Factors in Differentiating Adipocytes Disturbs the Expression of Autophagy Genes
We evaluated the role of CEBPB and PPARG in regulating the gene expression of key autophagy genes. We applied the gene differential expression analysis on two datasets of Cebpband Pparg-knockdown in adipocytes at different time point of the differentiation course ( Figure S8). The knockdown of either factor resulted in significant (log 2 fold-change < −75; FDR < 0.2) changes in the expression of lipogenic genes such as Acyl, Lpl and/or Fasn. These changes were pronounced in the late-adipocytes (day 5) with the knockdown of Pparg (Table 2 and Figure 9). The effect of the knockdown of both genes was also reflected by the enrichment of several lipid metabolism, transport and storage gene sets (Table S6). The disruption of the gene expression as a result of perturbing the adipogenic factor genes extended to other genes of interest.
For example, the expression of autophagy transcription factor gene Foxo1 was up-regulated (0.84 log 2 fold-change) by Cebpb-knockdown and down-regulated by Pparg-knockdown at day 2 (0.81 log 2 fold-change) (Table 2 and Figure 9). Conversely, Map1lc3 was down-regulated by Cebpb-knockdown (−0.54 log 2 fold-change) and up-regulated (0.93 log 2 fold-change) by Pparg-knockdown (Table 2 and Figure 9). In addition, several autophagy gene sets including those representing autophagosome assembly and maturation were enriched by the knockdown of the adipogenic factor genes (Table S6). Finally, the knockdown of Cebpb significantly reduced (<−0.65 log 2 fold-change) the expression of the other adipogenic transcription factor gene Pparg in the non-induced adipocyte and four hours after induction alike (Table 2 and Figure 9a). Removing Pparg on the other hand had only slight effect on the other adipogenic factor genes Cebpa and Cebpb and only on day 2 of the differentiation (Table 2 and Figure 9b).

Discussion
We found that, key autophagy markers such as Map1lc3b, Becn1 and Sqstm1 were regulated starting from the second day after induction of the adipocyte differentiation and reached the highest expression in the mature adipocytes. Some of these key genes were targeted directly by CEBPB and/or PPARG. Other key autophagy genes maybe regulated by the adipogenic transcription factors indirectly through other transcription factors such FOXO1, XBP1 and TFEB. Finally, the binding of the adipogenic factors on these genes is associated with the binding of other co-factors and histone modifications. Figure 10 depicts a summary diagram of the findings in this study. These findings were investigated under the knockdown of the adipogenic factor genes in differentiating adipocytes at the relevant time points.
Autophagy is important in the early stage of differentiation for lipid accumulation [2,10,11]. The knock-down of Atg7 impaired the differentiation process [2]. Our analysis suggests a progressive activation of autophagy genes starting at day 2 after induction as indicated by Map1lc3b, although still lower than its expression in the pre-adipocytes, which continues into the maturation stage ( Figure S3). Other studies suggested the activation of Atg4b through CEBPB [4] and the inhibition of SQSTM1 which enhances the ERK activation at day 2 [12,13]. Both events were observed in our analysis (Table 1). However, we did not observe a strong binding of CEBPB on any of the autophagy-related genes. Therefore, we investigated the binding of adipogenic factors on other key autophagy regulators.
Several studies suggested the direct regulation of autophagy by adipogenic transcription factors. The expression of Map1lc3b and Sqstm1 was directly controlled by CEBPB in hepatocytes [14]. CEBPG along with three other factors regulated the expression of CREB regulated transcription coactivator 2 (TORC2) in bovine adipocytes, which promotes autophagy [15,16]. Similarly, PPARG regulates autophagy through NEDD4 in HepG2 cells [5]. In agreement with the previous literature, we found that CEBPB binds around the promoter areas of Map1lc3b and Sqstm1 genes in addition to Becn1 although the binding intensity is more apparent after day 2 (Table S4 and Figure 4b,d). The latter also has binding sites for PPARG on which the occupancy of the transcription factor was the reverse of that of CEBPB over time (Figure 4b,c).
Other autophagy genes including autophagy-related genes that are under control of transcription factors TFEB, FOXO1 and/or XBP1 may be only indirectly regulated by adipogenic factors [17][18][19]. Both CEBPB and PPARG has binding sites on the genes encoding these three factors and appears to be actively regulated over time (Table S5 and Figure 5). Functionally, PPARG binding effect may be opposed to that of CEBPB on these genes as they have opposite binding patterns accompanied by up-regulation of their coding genes.  FOXO1 regulates adipocyte differentiation [7]. It inhibits the expression of Pparg therefore it should be deactivated early in the differentiation process to allow for the adipogenic factor activation [8]. Our analysis suggests that Foxo1 is turned on at some point in the differentiation course and continues to be expressed in the mature adipocytes (Table S5 and Figure 5a). This might be explained by the combined binding of the two adipogenic factors CEBPB and PPARG on the gene and their change of occupancy over the course of differentiation ( Figure 5). Here, the expression of PPARG plateaus and that of Foxo1 continues to rise.
CEBPB was reported to induce Pparg indirectly through CEBPA [20]. However, there seems to be an alternative mechanism to this induction as it happens even in the absence of the CEBP proteins [21,22]. Our analysis suggests multiple binding sites of CEBPB around the promoter of Pparg, this might explain the induction at lease in the case of absent CEBPA (Figure 6b,d). Together with the observed binding of PPARG on CEBPB forms a transcriptional loop that might enforce this induction mechanism (Figure 6b,c). Finally, CEBPB appears to bind to its own promoter to auto-regulate its own expression (Table S7 and Figure 6b,d). The binding of the two factors on the Cebpb promoter may explain the timely reduction in the expression of the gene as they reach their highest occupancy toward the end of the differentiation course.
A study reported a positive correlation between PPARG and the co-factor MED1 but not CBP in response to rosiglitazone in differentiating adipocytes [23]. Here, we investigated the possibility of this recruitment happening during the natural course of differentiation, whether it is shared by other transcription factors and co-factors. We found that the PPARG occupancy on autophagy genes correlates with the occupancy of MED1, RXRG, and EP300 ( Figure S6a). CEBPB occupancy only correlated with these factors in the pre-adipocytes ( Figure S6a). In both cases, the locations of the peaks had no effect on the correlation ( Figure S6b). The occupancy of the co-factors MED1 and RXRG correlated the most with PPARG while EP300 correlated with CEBPB across stages and genomic regions (Figures 7 and S7). However, the change in occupancy of PPARG and CEBPB on their own targets was positively correlated with the change of occupancy of MED1 and RXRG or MED1 and EP300 in early or late-differentiation respectively (Figure 8a,b). CEBPB and PPARG change in occupancy on their respective autophagy targets were mimicked by H3K27ac and H3K4me1 (Figure 8c,d). Both histone modifications were reported to mark active areas of the chromatin; the active genes in case of H3K27ac and active enhancers in the case of H3K4me1 [24,25]. This is consistent with the reported non-transcriptional role of PPARG which helps to reorganize the chromatin in the differentiating cells [26]. In other words, the areas of the chromatin which contain genes coding for autophagy proteins are activated as evidenced by the changes in the active chromatin markers (Table S4 and Figure 8b,c). These changes are associated with parallel changes in the adipogenic factor occupancy much like their non-transcriptional role in the spatial reorganization of adipogenic genes.
Although the knockdown of Cebpb or Pparg in the manner analyzed in this study may not reveal the functional effect on the phenotype, it provides confirmation of the mechanistic role of the transcription factors in regulating the expression of their own genes as well as of autophagy genes. Removing either transcription factors significantly reduced the expression of important lipogenic genes such as Acyl and Lpl (Table 2 and Figure 9). Similarly, the knockdown of Cebpb significantly reduced the expression of Pparg in non-induced and induced adipocytes (Table 2 and Figure 9a). As expected, the opposing regulation patterns of the adipogenic factors on key autophagy genes such as Foxo1 and Map1lc3b were reflected in their knockdown patterns ( Figure 9). However, we also expect some gene expression changes to be an indirect consequence of removing or relocating the factors to other regions. In addition, some important genes such as Becn1 and Tfeb didn't have probes corresponding to them in the microarrays dataset and we were not able to confirm their regulation under the Pparg-knockdown condition.
lThe 3T3-L1 cell line is a useful adipocyte model. It has many applications in obesity and insulin resistance research such as lipid synthesis, white vs. brown adipose tissue development, insulin-sensitizing drug action [27][28][29]. However, the findings based on this model should be considered to the extent that they reflect the development and the biology of the mature adipocytes. Therefore, the validation of such findings in primary human adipocytes or fat tissue isolates is important and is a natural extension to the current study. In addition, detailed experimentation may be required to study the conditions under which the regulatory links reported here result in qualitative changes in the gene expression and protein level of the targets. Finally, future studies are needed to confirm that the changes at the transcription level translate into autophagy activity.
Studying the reported findings under perturbation of the differentiation course would provide useful information. For example, it might be interesting to study the course of differentiation under autophagy inhibition or absence to quantify the changes in its gene expression profiles. Moreover, the conditional knockdown of CEBPB and/or PPARG in early and late-differentiation stages would provide details about the functional roles of the regulatory links reported in this study. Finally, the contributions of the direct and indirect regulation of the adipogenic transcription factors to autophagy genes can be determined by conditionally modifying or removing the intermediate autophagy transcription factors. Indeed, evaluating these findings under chemical and genetic modification of the differentiation course is a natural extension of this study and a goal for future investigation.
Using publicly available datasets from different sources represents a challenge to the current investigation. We addressed this challenge by manually curating the datasets, applying proper quality control measures and using proper quantitative analysis. Moreover, depending on the availability of the data, some of the phenotypes of interest might have a small number of samples, low time resolution or be missing altogether. In this study, we used quantitative models and reported all findings with an attached degree of uncertainty to reflect the degree of confidence in the data they were drawn from. Finally, we abstracted the time course of the differentiation into meaningful stages and key time points to address the missing data issues.

Pre-Adipocyte 3T3-L1 Differentiation Protocol
The data included in this study used the MDI (1-Methyl-3-isobutylxanthine, Dexamethasone, and Insulin) differentiation protocol with minimal variations [3]. The induction starts two days post-confluence and lasts for one or two weeks with changing the media at fixed intervals. At the end of it, the efficiency of the induction can be evaluated using adipogenic markers or lipid staining. The differentiation course was divided by time into three stages (non, 0 h and before; early, after 0 to 48 h; and late, after 48 to 260 h). The grouping criteria were previously suggested by others and devised from the data itself as it explains the largest amount of variance.

Data Collection, Pre-Processing, and Processing
The data collection strategy, pre-processing and processing of the datasets used in this study is documented in two Bioconductor packages curatedAdipoRNA and curatedAdipoChIP (submitted). Briefly, we surveyed the gene expression omnibus (GEO) and the sequence read archive (SRA) repositories for high-throughput sequencing data of MDI induced 3T3-L1 pre-adipocyte at different time points. The data were obtained from GEO/SRA in the form of raw reads. In total, 98 RNA-Seq and 207 transcription factor, co-factor and histone modification markers ChIP-Seq samples were included.
For RNA-Seq (n = 66, subset), the raw reads were aligned to mm10 mouse genome using HISAT2 [30]. featureCounts was used to count the reads aligned to the known genes [31]. For ChIP-Seq (n = 96, subset), the raw reads were aligned to mm10 using BOWTIE2, peaks were called and signal tracks were built using MACS2 and reads counted in a peak set of replicated peaks across samples were obtained using BEDTOOLS [32][33][34]. The peak set was annotated and peaks were assigned to the nearest gene using ChIPseeker [35]. FASTQC was used to assess the quality of the raw reads [36]. The phenotype data of the samples were manually curated to match the time point, stage of differentiation, the binding factors/markers and the ChIP antibodies (Tables S8 and S9).
Another dataset of CEBPB and PPARG ChIP-Seq samples in different tissue was used to test the binding specificity of the adipogenic transcription factors. The Cistrome Data Browser was searched for the adipocyte, fibroblast of macrophage (Biological Sources) and CEBPB/PPARG (Factors). The results were manually reviewed to select the ChIP-Seq samples (n = 11) of the two factors in these tissues (Table S10). The data were obtained in the form of processed signal tracks. Two datasets of 8 RNA-Seq and 18 microarrays samples of Cebpband Pparg-knockdown in differentiating adipocytes was obtained, processed and analyzed (Table S11). 3T3-L1 cells were transected with either shRNA against Cebpb or siRNA against Pparg. Transected cells were subsequently differentiated into mature adipocytes using the previously described protocol. Total RNA was harvested and the gene expression was profiled at zero/four hours or zero/two/five days post-induction using RNA-Seq or microarrays for Cebpband Pparg-knockdown cell, respectively.

Gene Expression Analysis
The gene counts were used to apply the differential expression analysis using DESeq2 [37]. The samples were divided into three groups corresponding to the stages none, early and late-differentiation. Three contrast groups were applied (early vs. non, late vs. non and late vs. early). For each gene, a log 2 fold-change and a false discovery rate (FDR) were calculated. Genes with absolute log 2 fold-change > 1 and FDR < 0.2 were considered significantly expressed and depending on the sign of the fold-change to be up-or down-regulated in one or more of the three contrasts. Using the gene counts a multidimensional scaling (MDS) analysis was applied to all, autophagy and adipogenic transcription factor genes using cmdscale (base R) [38]. The differential gene expression was applied to the Cebpb-knockdown RNA-Seq and Pparg-knockdown microarrays datasets by comparing the expression in knockdown (KD) vs. control at each time point using DESeq2 and limma [37,39].
The biological process gene ontology (GO) term autophagy (GO:0006914) was used to define the genes with known functions in autophagy (n = 158) [40]. The molecular function term DNA-binding transcription factor activity (GO:0003700) was used to define genes with known functions as transcription factors (n = 745). The intersection of the two terms were the six autophagy transcription factors genes.

Gene Set Enrichment Analysis
Several autophagy and lipid metabolism GO terms and their gene products were identified. The lists of differentially expressed genes were used to determine the over-represented gene sets in each comparison. For each gene set term, the ratio of the differentially expressed genes to the number of the genes in the term was compared to the total number of the differentially expressed genes to the total number of genes. A two-by-two table was constructed from the list and tested to calculate a p-value. The packages goseq and clusterProfiler were used to apply this analysis on the RNA-Seq and the microarrays datasets, respectively [43,44].

Gene Co-Expression Analysis
The gene counts were transformed using variance stabilization transformation (VST). The transformed counts were used to calculate the Pearson's correlation coefficient (PCC) as a measure of co-expression between each pair of genes in each condition. The coefficients were transformed into z-scores. The difference in z-score and the variance between every two conditions is used to calculate p-value. The z-scores from the original and permuted datasets are used to calculate empirical pand q-values for multiple testing adjustment. Pairs of genes were considered differentially co-expressed when FDR < 0.2. DGCA was used to apply this analysis [45].

Peak Binding Analysis
After removing the low-quality ChIP-Seq samples and the low count peaks in the peak set, the reads count in peaks were used to apply the differential peak binding analysis using DESeq2 [37]. The samples were divided into three groups based on the time point (non, 0 h and before; early, after 0 to 48 h; and late, after 48 to 260 h). Three contrast groups were applied (early vs. non, late vs. non and late vs. early). For each peak, a log 2 fold-change and an FDR were calculated. Peaks with absolute log 2 fold-change > 0.5 and FDR < 0.2 were considered significantly expressed and depending on the sign of the fold-change to be up-or down-regulated in one or more of the three contrasts.

Occupancy and Affinity Analyses
The reads count in peaks were aggregated by gene to find the total occupancy of the factors and histone modification markers. The occupancy of the factors and histone markers were grouped by genomic annotations (Promoter, 3' Untranslated Region (UTR), 5' UTR or other). The occupancy correlation of each sample was calculated using the PCC of all pairs reads count in peaks, total occupancy or occupancy grouped by genomic annotation. Promoters were defined as ±3 kb around the transcription start site (TSS) of each transcript. The signal tracks were used to visualize the peaks in the promoter regions of selected genes as an enrichment score relative to known controls.

Source Code and Reproducibility
The input, output and the tools used in each step of the study workflow is depicted in Figure 1. The software environment where the full analysis and manuscript production were done is available as a docker container for reproducibility (https://hub.docker.com/r/bcmslab/autoreg). The analysis was conducted in R (3.5) using Bioconductor (3.6) and other R packages referred to in the relevant subsections [46,58]. The source code for conducting this analysis is maintained on (https://github. com/BCMSLab/autoreg). The code for reproducing the figures and tables presented in this manuscript is available at (https://github.com/BCMSLab/auto_adipo_diff).

Conclusions
Autophagy genes are regulated as part of the differentiation course of the adipocytes. This regulation is driven by adipogenic transcription factors such as CEBPB and PPARG. The adipogenic factors target key autophagy genes such as Becn1, Map1lc3b and Sqstm1. Moreover, other autophagy genes are regulated indirectly through autophagy transcription factors Foxo1, Tfeb and Xbp1. Co-factors such as MED1 actively contribute to the transcription factors binding on autophagy genes. Other co-factors such as RXRG and EP300 associate specifically with PPARG or CEBPB respectively. H3K4me1 and H3K27ac markers also associate with the adipogenic factors binding sites indicating non-transcriptional roles such as organizing the chromatin regions containing autophagy genes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4409/8/11/1321/ s1, Table S1. Percent of gene expression variance explained by the stage of differentiation; Table S2. Differential expression of gene markers; Table S3. Autophagy and lipid metabolism gene ontology terms enrichment in differentiating adipocytes; Table S4. Significant peaks of adipogenic factors on autophagy genes; Table S5. Significant peaks of adipogenic factors on autophagy transcription factor genes; Table S6. Autophagy and lipid metabolism gene ontology terms enrichment in Cebpband Pparg-knockdown differentiating adipocytes; Table S7. Peaks of adipogenic factors on adipogenic transcription factor genes; Table S8. Datasets of RNA and ChIP Seq; Table S9. ChIP antibodies for transcription factors, co-factors and histone markers; Table S10. Datasets of ChIP Seq in different tissues; Table S11. Datasets of RNA Seq and microarrays of adipocytes with Cebpb or Pparg-knockdown; Figure S1. Differential gene expression in early and late differentiating adipocytes; Figure S2. Differential peak binding of adipogenic transcription factors on autophagy genes; Figure S3. Gene markers of differentiation, lipogenesis and autophagy; Figure S4. Co-expression of adipogenic transcription factors genes with autophagy genes; Figure S5. Tissue specific binding of adipogenic transcription factors; Figure S6. Adipogenic transcription factors correlations with co-factors and histone markers on autophagy genes; Figure S7. Correlation in occupancy between adipogenic transcription factors and histone markers on autophagy genes at different genomic locations; Figure S8. Differential gene expression in Cebpb or Pparg-knockdown differentiating adipocytes.

TSS
Transcription Start Site UTR Un-translated Region