Next Article in Journal
Perturbation of Methionine/S-adenosylmethionine Metabolism as a Novel Vulnerability in MLL Rearranged Leukemia
Next Article in Special Issue
Autophagy Attenuation Hampers Progesterone Synthesis during the Development of Pregnant Corpus Luteum
Previous Article in Journal
Wnt Glycation Inhibits Canonical Signaling
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cells
Volume 8
Issue 11
10.3390/cells8111321
cells-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Mahmoud Ahmed
,
Trang Huyen Lai
,
Jin Seok Hwang
,
Sahib Zada
,
Trang Minh Pham
and
Deok Ryong Kim
*
Department of Biochemistry and Convergence Medical Sciences and Institute of Health Sciences, Gyeongsang National University School of Medicine, Jinju 527-27, Korea
*
Author to whom correspondence should be addressed.
Cells 2019, 8(11), 1321; https://doi.org/10.3390/cells8111321
Submission received: 4 October 2019 / Revised: 20 October 2019 / Accepted: 22 October 2019 / Published: 25 October 2019
(This article belongs to the Special Issue Lipids and Lipid Metabolism in Autophagy)
Browse Figures
Review Reports Versions Notes

Abstract

:
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.

1. 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].
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.

2. Results

2.1. 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).

2.2. 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).

2.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).

2.4. 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).

2.5. 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).

2.6. 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 (Figure 7 and Figure S7). On the other hand, CEBPB is closest to the co-factor E1A binding protein p300 (EP300) at the same stage and location (Figure 7 and Figure 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).

2.7. 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 Cebpb- and 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).

3. 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 (Figure 7 and Figure 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.

4. Materials and Methods

4.1. 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.

4.2. 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 Cebpb- and 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 Cebpb- and Pparg-knockdown cell, respectively.

4.3. 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. The GO annotations were accessed using org.Mm.eg.db [41]. Other genomic annotation such as gene coordinates and identifiers were accessed using TxDb.Mmusculus.UCSC.mm10.knownGene [42].

4.4. 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].

4.5. 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 p- and 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].

4.6. 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.

4.7. 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.

4.8. Data Management, Transformation and Visualization

The processed data were stored mainly in the ExpressionSet, SummarizedExperiment and bigwig formats [46,47]. GEOquery was used to obtain the microarrays data [48]. collapseRows (WGCNA), tidyverse, GenomicRanges and rtracklayer were used for data transformation [49,50,51,52]. The ggplot2, GGally, xtable, ComplexHeatmap and EnrichedHeatmap were used for data visualization [53,54,55,56,57].

4.9. 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).

5. 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 https://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 Cebpb- and 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.

Author Contributions

M.A. collected, processed and analyzed the data and prepared the manuscript. T.H.L., J.S.H., S.Z. and T.M.P. contributed to preparing the manuscript. D.R.K. Supervised the study and contributed to writing and editing the manuscript.

Funding

This study was supported by the Basic Research Program of the National Research Foundation of Korea (NRF) and founded by the Ministry of Education, Science and Technology (2018R1D1A1B07043715) and the Ministry of Science, ICT and Future Planning (2015R1A5A2008833).

Acknowledgments

We thank the lab members for the critical discussion of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
FDRFalse Discovery Rate
GEOGene Expression Omnibus
GOGene Ontology
MDI1-Methyl-3-isobutylxanthine, Dexamethasone and Insulin
MDSMultidimensional Scaling
PCCPearson’s Correlation Coefficient
SRASequence Read Archive
TSSTranscription Start Site
UTRUn-translated Region

References

  1. Dong, H.; Czaja, M.J. Regulation of lipid droplets by autophagy. Trends Endocrinol. Metab. TEM 2011, 22, 234–240. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Singh, R.; Xiang, Y.; Wang, Y.; Baikati, K.; Cuervo, A.M.; Luu, Y.K.; Tang, Y.; Pessin, J.E.; Schwartz, G.J.; Czaja, M.J. Autophagy regulates adipose mass and differentiation in mice. J. Clin. Investig. 2009, 119, 3329–3339. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Green, H.; Kehinde, O. An established preadipose cell line and its differentiation in culture II. Factors affecting the adipose conversion. Cell 1975, 5, 19–27. [Google Scholar] [CrossRef]
  4. Guo, L.; Huang, J.X.; Liu, Y.Y.; Li, X.; Zhou, S.R.; Qian, S.W.; Liu, Y.Y.; Zhu, H.; Huang, H.Y.; Dang, Y.J.; et al. Transactivation of Atg4b by C/EBPβ promotes autophagy to facilitate adipogenesis. Mol. Cell. Biol. 2013, 33, 3180–3190. [Google Scholar] [CrossRef]
  5. Liu, J.; Yao, Q.; Xiao, L.; Ma, W.; Li, F.; Lai, B.; Wang, N. PPAR γ induces NEDD 4 gene expression to promote autophagy and insulin action. FEBS J. 2019. [Google Scholar] [CrossRef]
  6. Zhao, Y.; Yang, J.; Liao, W.; Liu, X.; Zhang, H.; Wang, S.; Wang, D.; Feng, J.; Yu, L.; Zhu, W.G. Cytosolic FoxO1 is essential for the induction of autophagy and tumour suppressor activity. Nat. Cell Biol. 2010, 12, 665–675. [Google Scholar] [CrossRef]
  7. Nakae, J.; Kitamura, T.; Kitamura, Y.; Biggs, W.H.; Arden, K.C.; Accili, D. The forkhead transcription factor Foxo1 regulates adipocyte differentiation. Dev. Cell 2003, 4, 119–129. [Google Scholar] [CrossRef]
  8. Armoni, M.; Harel, C.; Karni, S.; Chen, H.; Bar-Yoseph, F.; Ver, M.R.; Quon, M.J.; Karnieli, E. FOXO1 Represses Peroxisome Proliferator-activated Receptor-γ1 and -γ2 Gene Promoters in Primary Adipocytes. J. Biol. Chem. 2006, 281, 19881–19891. [Google Scholar] [CrossRef]
  9. Ahmed, M.; Quoc Nguyen, H.; Seok Hwang, J.; Zada, S.; Huyen Lai, T.; Soo Kang, S.; Ryong Kim, D. Systematic characterization of autophagy-related genes during the adipocyte differentiation using public-access data. Oncotarget 2018, 9. [Google Scholar] [CrossRef]
  10. Baerga, R.; Zhang, Y.; Chen, P.H.; Goldman, S.; Jin, S. Targeted deletion of autophagy-related 5 (atg5) impairs adipogenesis in a cellular model and in mice. Autophagy 2009, 5, 1118–1130. [Google Scholar] [CrossRef]
  11. Zhang, Y.; Goldman, S.; Baerga, R.; Zhao, Y.; Komatsu, M.; Jin, S. Adipose-specific deletion of autophagy- related gene 7 (atg7) in mice reveals a role in adipogenesis. Proc. Natl. Acad. Sci. USA 2009, 106, 19860–19865. [Google Scholar] [CrossRef] [PubMed]
  12. Rodriguez, A.; Durán, A.; Selloum, M.; Champy, M.F.; Diez-Guerra, F.J.; Flores, J.M.; Serrano, M.; Auwerx, J.; Diaz-Meco, M.T.; Moscat, J. Mature-onset obesity and insulin resistance in mice deficient in the signaling adapter p62. Cell Metab. 2006, 3, 211–222. [Google Scholar] [CrossRef] [PubMed]
  13. Lee, S.J.; Pfluger, P.T.; Kim, J.Y.; Nogueiras, R.; Duran, A.; Pagès, G.; Pouysségur, J.; Tschöp, M.H.; Diaz-Meco, M.T.; Moscat, J. A functional role for the p62-ERK1 axis in the control of energy homeostasis and adipogenesis. EMBO Rep. 2010, 11, 226–232. [Google Scholar] [CrossRef] [PubMed]
  14. Ma, D.; Panda, S.; Lin, J.D. Temporal orchestration of circadian autophagy rhythm by C/EBPβ. EMBO J. 2011, 30, 4642–4651. [Google Scholar] [CrossRef]
  15. Khan, R.; Raza, S.H.A.; Junjvlieke, Z.; Wang, X.; Garcia, M.; Elnour, I.E.; Wang, H.; Linsen, Z. Function and Transcriptional Regulation of Bovine TORC2 Gene in Adipocytes: Roles of C/EBP, XBP1, INSM1 and ZNF263. Int. J. Mol. Sci. 2019, 20, 4338. [Google Scholar] [CrossRef] [PubMed]
  16. Vlahakis, A.; Powers, T. A role for TOR complex 2 signaling in promoting autophagy. Autophagy 2014. [Google Scholar] [CrossRef]
  17. Settembre, C.; Di Malta, C.; Polito, V.A.; Garcia Arencibia, M.; Vetrini, F.; Erdin, S.S.U.; Erdin, S.S.U.; Huynh, T.; Medina, D.; Colella, P.; et al. TFEB links autophagy to lysosomal biogenesis. Science 2011, 332, 1429–1433. [Google Scholar] [CrossRef]
  18. Fiorentino, L.; Cavalera, M.; Menini, S.; Marchetti, V.; Mavilio, M.; Fabrizi, M.; Conserva, F.; Casagrande, V.; Menghini, R.; Pontrelli, P.; et al. Loss of TIMP3 underlies diabetic nephropathy via FoxO1/STAT1 interplay. EMBO Mol. Med. 2013, 5, 441–455. [Google Scholar] [CrossRef]
  19. Margariti, A.; Li, H.; Chen, T.; Martin, D.; Vizcay-Barrena, G.; Alam, S.; Karamariti, E.; Xiao, Q.; Zampetaki, A.; Zhang, Z.; et al. XBP1 mRNA splicing triggers an autophagic response in endothelial cells through BECLIN-1 transcriptional activation. J. Biol. Chem. 2013, 288, 859–872. [Google Scholar] [CrossRef]
  20. Wu, Z.; Bucher, N.L.; Farmer, S.R. Induction of peroxisome proliferator-activated receptor gamma during the conversion of 3T3 fibroblasts into adipocytes is mediated by {C}/{EBPbeta}, {C}/{EBPdelta}, and glucocorticoids. Mol. Cell. Biol. 1996, 16, 4128–4136. [Google Scholar] [CrossRef]
  21. Wu, Z.; Rosen, E.D.; Brun, R.; Hauser, S.; Adelmant, G.; Troy, A.E.; McKeon, C.; Darlington, G.J.; Spiegelman, B.M. Cross-regulation of C/EBP alpha and PPAR gamma controls the transcriptional pathway of adipogenesis and insulin sensitivity. Mol. Cell 1999, 3, 151–158. [Google Scholar] [CrossRef]
  22. Tanaka, T.; Yoshida, N.; Kishimoto, T.; Akira, S. Defective adipocyte differentiation in mice lacking the C/EBPβ and/or C/EBPδ gene. EMBO J. 1997, 16, 7432–7443. [Google Scholar] [CrossRef] [PubMed]
  23. Haakonsson, A.K.; Stahl Madsen, M.; Nielsen, R.; Sandelin, A.; Mandrup, S. Acute Genome-Wide Effects of Rosiglitazone on PPARγ Transcriptional Networks in Adipocytes. Mol. Endocrinol. 2013, 27, 1536–1549. [Google Scholar] [CrossRef] [PubMed]
  24. Creyghton, M.P.; Cheng, A.W.; Welstead, G.G.; Kooistra, T.; Carey, B.W.; Steine, E.J.; Hanna, J.; Lodato, M.A.; Frampton, G.M.; Sharp, P.A.; et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. USA 2010, 107, 21931–21936. [Google Scholar] [CrossRef] [Green Version]
  25. Rada-Iglesias, A. Is H3K4me1 at enhancers correlative or causative? Nat. Genet. 2018, 50, 4–5. [Google Scholar] [CrossRef]
  26. Sarusi Portuguez, A.; Schwartz, M.; Siersbaek, R.; Nielsen, R.; Sung, M.H.; Mandrup, S.; Kaplan, T.; Hakim, O. Hierarchical role for transcription factors and chromatin structure in genome organization along adipogenesis. FEBS J. 2017, 284, 3230–3244. [Google Scholar] [CrossRef]
  27. Glenn, K.C.; Shieh, J.J.; Laird, D.M. Characterization of 3T3-L1 storage lipid metabolism: Effect of somatotropin and insulin on specific pathways. Endocrinology 1992. [Google Scholar] [CrossRef]
  28. Sarjeant, K.; Stephens, J.M. Adipogenesis. Cold Spring Harb. Perspect. Biol. 2012, 4, a008417. [Google Scholar] [CrossRef]
  29. Thomson, M.J.; Williams, M.G.; Frost, S.C. Development of insulin resistance in 3T3-L1 adipocytes. J. Biol. Chem. 1997. [Google Scholar] [CrossRef]
  30. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef]
  31. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930. [Google Scholar] [CrossRef] [PubMed]
  32. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10, R25. [Google Scholar] [CrossRef] [PubMed]
  33. Zhang, Y.; Liu, T.; Meyer, C.A.; Eeckhoute, J.; Johnson, D.S.; Bernstein, B.E.; Nusbaum, C.; Myers, R.M.; Brown, M.; Li, W.; et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008, 9, R137. [Google Scholar] [CrossRef] [PubMed]
  34. Quinlan, A.R.; Hall, I.M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26, 841–842. [Google Scholar] [CrossRef] [PubMed]
  35. Yu, G.; Wang, L.G.; He, Q.Y. ChIP seeker: An R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 2015, 31, 2382–2383. [Google Scholar] [CrossRef]
  36. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. Java Package. 2018. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed on 4 October 2019).
  37. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  38. Cox, T.F.; Cox, M. Multidimensional Scaling, 2nd ed.; Chapman and Hall: New York, NY, USA, 2000; p. 328. [Google Scholar] [CrossRef]
  39. Smyth, G.K. Limma: Linear Models for Microarray Data. In Bioinformatics and Computational Biology Solutions Using R and Bioconductor; Gentleman, R., Carey, V.J., Huber, W., Irizarry, R.A., Eds.; Springer: New York, NY, USA, 2005; Chapter 23; pp. 397–420. [Google Scholar]
  40. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene Ontology: Tool for the unification of biology. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef]
  41. Carlson, M. org.Mm.eg.db: Genome Wide Annotation for Mouse. R Package Version 3.8.2. 2019. Available online: https://doi.org/doi:10.18129/B9.bioc.org.Mm.eg.db (accessed on 4 October 2019).
  42. Team, B.C.; Maintainer, B.P. TxDb.Mmusculus.UCSC.mm10.knownGene: Annotation Package for TxDb Object(s). R Package Version 3.4.7. 2019. Available online: https://doi.org/doi:10.18129/B9.bioc.TxDb.Mmusculus.UCSC.mm10.knownGene (accessed on 4 October 2019).
  43. Young, M.D.; Wakefield, M.J.; Smyth, G.K.; Oshlack, A. Gene ontology analysis for RNA-seq: Accounting for selection bias. Genome Biol. 2010, 11, R14. [Google Scholar] [CrossRef]
  44. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef]
  45. McKenzie, A.T.; Katsyv, I.; Song, W.M.; Wang, M.; Zhang, B. DGCA: A comprehensive R package for Differential Gene Correlation Analysis. BMC Syst. Biol. 2016, 10, 1–25. [Google Scholar] [CrossRef] [Green Version]
  46. Huber, W.; Carey, V.J.; Gentleman, R.; Anders, S.; Carlson, M.; Carvalho, B.S.; Bravo, H.C.; Davis, S.; Gatto, L.; Girke, T.; et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods 2015, 12, 115–121. [Google Scholar] [CrossRef] [PubMed]
  47. Morgan, M.; Obenchain, V.; Hester, J.; Pagès, H. SummarizedExperiment: SummarizedExperiment Container. R Package Version 1.14.1. 2019. Available online: https://doi.org/doi:10.18129/B9.bioc.SummarizedExperiment (accessed on 4 October 2019).
  48. Sean, D.; Meltzer, P.S. GEOquery: A bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007, 23, 1846–1847. [Google Scholar] [CrossRef]
  49. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef] [PubMed]
  50. Wickham, H. tidyverse: Easily Install and Load the “Tidyverse”. R Package Version 1.2.1. 2017. Available online: https://cran.r-project.org/package=tidyverse (accessed on 4 October 2019).
  51. Lawrence, M.; Huber, W.; Pagès, H.; Aboyoun, P.; Carlson, M.; Gentleman, R.; Morgan, M.T.; Carey, V.J. Software for Computing and Annotating Genomic Ranges. PLoS Comput. Biol. 2013, 9, e1003118. [Google Scholar] [CrossRef]
  52. Lawrence, M.; Gentleman, R.; Carey, V. rtracklayer: An R package for interfacing with genome browsers. Bioinformatics 2009, 25, 1841–1842. [Google Scholar] [CrossRef]
  53. Wickham, H. Ggplot2; Springer: New York, NY, USA, 2011. [Google Scholar] [CrossRef]
  54. Schloerke, B.; Crowley, J.; Cook, D.; Briatte, F.; Marbach, M.; Thoen, E.; Elberg, A.; Larmarange, J. GGally: Extension to “ggplot2”. R Package Version 1.4.0. 2018. Available online: https://cran.r-project.org/package=GGally (accessed on 4 October 2019).
  55. Dahl, D.B.; Scott, D.; Roosen, C.; Magnusson, A.; Swinton, J. xtable: Export Tables to LaTeX or HTML. R Package Version 1.8-4. 2018. Available online: https://cran.r-project.org/package=xtable (accessed on 4 October 2019).
  56. Gu, Z.; Eils, R.; Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016, 32, 2847–2849. [Google Scholar] [CrossRef] [Green Version]
  57. Gu, Z.; Eils, R.; Schlesner, M.; Ishaque, N. EnrichedHeatmap: An R/Bioconductor package for comprehensive visualization of genomic signal associations. BMC Genom. 2018, 19, 234. [Google Scholar] [CrossRef]
  58. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017. [Google Scholar]
Cells 08 01321 g001 550
Figure 1. Data pre-processing, processing and analyses’ workflow of RNA-Seq and ChIP-Seq datasets.Pre-processing, processing and analyses’ workflow
Figure 1. Data pre-processing, processing and analyses’ workflow of RNA-Seq and ChIP-Seq datasets.Pre-processing, processing and analyses’ workflow
Cells 08 01321 g001
Cells 08 01321 g002 550
Figure 2. Multidimensional scaling analysis of gene expression. The transformed read counts of 15,786 genes coding, 158 autophagy genes, six genes encoding autophagy transcription factors (Foxo1, Irf8, Tfeb, Trp53, Xbp1 and Zkscan3) and five genes encoding adipogenic transcription factors (Cebpb, Ep300, Med1, Pparg and Rxrg) were used to perform multidimensional scaling (MDS). Counts were extracted from RNA-Seq samples (n = 66) of MDI-induced 3T3-L1 cells in three differentiation stages; 18 non (red), 25 early (green) and 23 late (blue) differentiating samples.
Figure 2. Multidimensional scaling analysis of gene expression. The transformed read counts of 15,786 genes coding, 158 autophagy genes, six genes encoding autophagy transcription factors (Foxo1, Irf8, Tfeb, Trp53, Xbp1 and Zkscan3) and five genes encoding adipogenic transcription factors (Cebpb, Ep300, Med1, Pparg and Rxrg) were used to perform multidimensional scaling (MDS). Counts were extracted from RNA-Seq samples (n = 66) of MDI-induced 3T3-L1 cells in three differentiation stages; 18 non (red), 25 early (green) and 23 late (blue) differentiating samples.
Cells 08 01321 g002
Cells 08 01321 g003 550
Figure 3. Co-expression of autophagy products with adipogenic transcription factors. The read counts of five genes coding for adipogenic transcription factors/co-factors and 158 autophagy genes were used to calculate the co-expression values of each of the factors with autophagy genes. Pearson’s correlation co-efficient were calculated from RNA-Seq samples (n = 66) of MDI-induced 3T3-L1 cells in three differentiation stages; non, early and late differentiating cell. Genes were stratified by their differential expression into down- (red), none- (green) or up-regulated (blue).
Figure 3. Co-expression of autophagy products with adipogenic transcription factors. The read counts of five genes coding for adipogenic transcription factors/co-factors and 158 autophagy genes were used to calculate the co-expression values of each of the factors with autophagy genes. Pearson’s correlation co-efficient were calculated from RNA-Seq samples (n = 66) of MDI-induced 3T3-L1 cells in three differentiation stages; non, early and late differentiating cell. Genes were stratified by their differential expression into down- (red), none- (green) or up-regulated (blue).
Cells 08 01321 g003
Cells 08 01321 g004 550
Figure 4. Key autophagy gene targets of the adipogenic transcription factors. (a) The reads count of three key autophagy genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of key autophagy genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples (n = 9 and 13) was used to represent peak enrichment. The Peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Figure 4. Key autophagy gene targets of the adipogenic transcription factors. (a) The reads count of three key autophagy genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of key autophagy genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples (n = 9 and 13) was used to represent peak enrichment. The Peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Cells 08 01321 g004
Cells 08 01321 g005 550
Figure 5. Autophagy transcription factor gene targets of the adipogenic transcription factors. (a) The reads count of three autophagy transcription factor coding genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of three autophagy transcription factor coding genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples was used to represent peak enrichment. The Peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Figure 5. Autophagy transcription factor gene targets of the adipogenic transcription factors. (a) The reads count of three autophagy transcription factor coding genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of three autophagy transcription factor coding genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples was used to represent peak enrichment. The Peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Cells 08 01321 g005
Cells 08 01321 g006 550
Figure 6. Auto-regulation and feedback loops of the adipogenic transcription factors. (a) The reads count of two adipogenic transcription factor coding genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of two adipogenic transcription factor coding genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples was used to represent peak enrichment. The peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Figure 6. Auto-regulation and feedback loops of the adipogenic transcription factors. (a) The reads count of two adipogenic transcription factor coding genes in RNA-Seq samples (n = 66) were used to represent the gene expression. The reads count in peaks of POLR2A in genomic region of two adipogenic transcription factor coding genes from ChIP-Seq samples (n = 12) was used to represent the transcription activity. (b) The reads count in peaks of CEBPB and PPARG ChIP-Seq samples was used to represent peak enrichment. The peaks from the (c) PPARG and (d) CEBPB ChIP-Seq samples were visualized by the time point and position around the transcription start site ( ± 3 kb). Cidec and klf5 were used as positive control genes.
Cells 08 01321 g006
Cells 08 01321 g007 550
Figure 7. 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.
Figure 7. 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.
Cells 08 01321 g007
Cells 08 01321 g008 550
Figure 8. Change in occupancy of transcription factors and co-factors and histone marker tags. The total count of reads in peaks of autophagy genes (n = 158) of ChIP-Seq samples (n = 75) of transcription factors, co-factors and histone modification was used calculate the occupancy. The change in occupancy of (a,c) PPARG and (b,d) CEBPB in early or late differentiation vs. the non-differentiated samples were plotted against the change in occupancy of co-factors (a,b) or histone modification tags (c,d). Points represent the fold-change in occupancy between stages of differentiation.
Figure 8. Change in occupancy of transcription factors and co-factors and histone marker tags. The total count of reads in peaks of autophagy genes (n = 158) of ChIP-Seq samples (n = 75) of transcription factors, co-factors and histone modification was used calculate the occupancy. The change in occupancy of (a,c) PPARG and (b,d) CEBPB in early or late differentiation vs. the non-differentiated samples were plotted against the change in occupancy of co-factors (a,b) or histone modification tags (c,d). Points represent the fold-change in occupancy between stages of differentiation.
Cells 08 01321 g008
Cells 08 01321 g009 550
Figure 9. Expression of key adipogenic, autophagy and lipogenic genes in Cebpb or Pparg-knockdown adipocytes. (a) Reads count of selected genes from RNA-Seq samples (n = 8) of Cebpb-knockdown differentiating adipocytes at 0 and 4 h were scaled and shown as color values (low, green; high, blue). (b) Probe intensity of selected genes from microarray samples (n = 18) of Pparg-knockdown differentiating adipocytes at 0, 2 and 5 days were scaled and shown as color values (low, green; high, blue). Genes are indicated by symbols in rows and grouped in four categories (Adipogenic TF, Autophagy Gene, Autophagy TF or Lipogenesis). Sample group (knockdown vs. control) is indicated in rows and grouped by time points (0 or 4 h; and 0, 2 or 5 days).
Figure 9. Expression of key adipogenic, autophagy and lipogenic genes in Cebpb or Pparg-knockdown adipocytes. (a) Reads count of selected genes from RNA-Seq samples (n = 8) of Cebpb-knockdown differentiating adipocytes at 0 and 4 h were scaled and shown as color values (low, green; high, blue). (b) Probe intensity of selected genes from microarray samples (n = 18) of Pparg-knockdown differentiating adipocytes at 0, 2 and 5 days were scaled and shown as color values (low, green; high, blue). Genes are indicated by symbols in rows and grouped in four categories (Adipogenic TF, Autophagy Gene, Autophagy TF or Lipogenesis). Sample group (knockdown vs. control) is indicated in rows and grouped by time points (0 or 4 h; and 0, 2 or 5 days).
Cells 08 01321 g009
Cells 08 01321 g010 550
Figure 10. Direct and indirect regulation of autophagy genes by adipogenic transcription factors. A constructed model of the transcriptional activity of two adipogenic factors (CEBPB and PPARG) on their corresponding genes, key autophagy genes (Map1lc3, Becn1 and Sqstm1) and autophagy transcription factors (Xbp1, Foxo1 and Tfeb). The direction of regulation (←, induction; or ⊢, repression) is inferred based on the availability of the factor in each stage and the co-expression with the target gene. The change in occupancy and the binding intensity is indicated by the edge type (─, high; or ─, low).
Figure 10. Direct and indirect regulation of autophagy genes by adipogenic transcription factors. A constructed model of the transcriptional activity of two adipogenic factors (CEBPB and PPARG) on their corresponding genes, key autophagy genes (Map1lc3, Becn1 and Sqstm1) and autophagy transcription factors (Xbp1, Foxo1 and Tfeb). The direction of regulation (←, induction; or ⊢, repression) is inferred based on the availability of the factor in each stage and the co-expression with the target gene. The change in occupancy and the binding intensity is indicated by the edge type (─, high; or ─, low).
Cells 08 01321 g010
Table
Table 1. Significant differentially expressed genes of adipogenic and autophagy transcription factors. FC, fold-change; SE, standard error.
Table 1. Significant differentially expressed genes of adipogenic and autophagy transcription factors. FC, fold-change; SE, standard error.
CategoryGeneEarly vs. NonLate vs. NonLate vs. Early
FCSEFCSEFCSE
Adipogenic TFCebpb1.50.19 −1.30.18
Med10.30.10.280.1
Pparg1.550.252.760.251.210.23
Rxrg 7.890.726.890.63
Autophagy TFFoxo11.10.222.030.220.920.2
Tfeb 1.660.261.650.24
Trp53−0.520.18−1.250.19−0.730.17
Xbp11.440.171.870.170.440.16
Zkscan3 0.590.120.410.11
Autophagy GeneAtg4b0.380.07 −0.460.06
Becn10.310.080.440.090.130.08
Map1lc3a−0.50.14 0.620.13
Map1lc3b−0.720.17−0.270.170.450.15
Sqstm1−1.020.15 0.910.14
Ulk1 0.740.150.550.14
Table
Table 2. Significant differentially expressed genes in Cebpb or Pparg knockdown adipocytes. KD, knockdown; h, hour; d, day; FC, fold-change; SE, standard error.
Table 2. Significant differentially expressed genes in Cebpb or Pparg knockdown adipocytes. KD, knockdown; h, hour; d, day; FC, fold-change; SE, standard error.
CategoryGeneCebpb KD vs. ControlPparg KD vs. Control
0 h4 h0 d2 d5 d
FCSEFCSEFCSEFCSEFCSE
Adipogenic TFCebpb−2.610.17−1.940.17 −0.280.11
Pparg−0.650.17−1.190.27−2.270.18−1.920.19−2.530.1
Cebpa 0.290.11
Autophagy GeneMap1lc3b−0.540.12−0.670.19 0.930.3
Sqstm1−0.580.11
Autophagy TFFoxo10.410.190.840.2 −0.810.26
Xbp1−0.30.15 −0.880.09−0.620.08
LipogenesisAcly−0.750.11 −0.740.14
Lpl−2.540.37−2.250.19−2.180.21−1.720.08−1.360.1
Fasn −0.910.31

Share and Cite

MDPI and ACS Style

Ahmed, M.; Lai, T.H.; Hwang, J.S.; Zada, S.; Pham, T.M.; Kim, D.R. 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. Cells 2019, 8, 1321. https://doi.org/10.3390/cells8111321

AMA Style

Ahmed M, Lai TH, Hwang JS, Zada S, Pham TM, Kim DR. 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. Cells. 2019; 8(11):1321. https://doi.org/10.3390/cells8111321

Chicago/Turabian Style

Ahmed, Mahmoud, Trang Huyen Lai, Jin Seok Hwang, Sahib Zada, Trang Minh Pham, and Deok Ryong Kim. 2019. "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" Cells 8, no. 11: 1321. https://doi.org/10.3390/cells8111321

APA Style

Ahmed, M., Lai, T. H., Hwang, J. S., Zada, S., Pham, T. M., & Kim, D. R. (2019). 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. Cells, 8(11), 1321. https://doi.org/10.3390/cells8111321

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop