Systematic characterization of autophagy-related genes during the adipocyte differentiation using public-access data

Autophagy contributes to reorganizing intracellular components and forming fat droplets during the adipocyte differentiation. Here, we systematically describe the role of autophagy-related genes and gene sets during the differentiation of adipocytes. We used a public dataset from the European Nucleotide Archive from an RNA-seq experiment in which 3T3-L1 cells were induced by a differentiation induction medium, total RNA was extracted and sequenced at four different time points. Raw reads were aligned to the UCSC mouse reference genome (mm10) using HISAT2, and aligned reads were summarized at the gene or exon level using HTSeq. DESeq2 and DEXSeq were used to model the gene and exon counts and test for differential expression and relative exon usage, respectively. After applying the appropriate transformation, gene counts were used to perform the gene set and pathway enrichment analysis. Data were obtained, processed and annotated using R and Bioconductor. Several autophagy-related genes and autophagy gene sets, as defined in the Gene Ontology, were actively regulated during the course of the adipocyte differentiation. We further characterized these gene sets by clustering their members to a few distinct temporal profiles. Other potential functionally related genes were identified using a machine learning procedure. In summary, we characterized the autophagy gene sets and their members to biologically meaningful groups and elected a number of genes to be functionally related based on their expression patterns, suggesting that autophagy plays a critical role in removal of some intracellular components and supply of energy sources for lipid biogenesis during adipogenesis.


INTRODUCTION
3T3-L1 pre-adipocyte is a mouse fibroblast that has the potential to differentiate into a mature adipocyte when treated with the MDI (160 nM insulin, 250 nM dexamethasone, and 0.5 mM 1-methyl-3-isobutylxanthine) differentiation induction medium [1]. Therefore, MDIinduced 3T3-L1 cells can be used as a platform for studying the adipocyte differentiation and further understanding the development of obesity and insulin resistance. Autophagy, a lysosome-dependent degradation process, is involved in many cellular mechanisms including stress response, cell growth and death, and cell differentiation. Particularly, it contributes to reorganizing intracellular components and forming fat droplets during the adipocyte differentiation [2]. Indeed, inhibition of autophagy during the adipogenic induction of 3T3-L1 cells results in the prevention of formation of lipid droplets [3].
High through-put technologies such as microarrays and RNA-seq facilitate the study of the cell biology by enabling the investigation of a large number of genetic elements simultaneously compared to the conventional Research Paper www.oncotarget.com laboratory techniques. The availability of public-access data in repositories such as the Sequence Read Archive (SRA) and the European Nucleotide Archive (ENA) enables biologists to reexamine previously published data for their own purposes of discovery and hypothesis generation [4,5].
Here, we use a public-access RNA-seq dataset of MDI-induced 3T3-L1 pre-adipocyte to study the role of autophagy in the adipocyte differentiation. First, we examine the changes in gene expression and relative exon usage over time of autophagy-related genes-members of the Gene Ontology term autophagy (GO:0006914). Then, we perform a gene set enrichment analysis and visualize the autophagy pathways to exploit the potential changes and modification in differentiating adipocytes compared to pre-adipocytes. The overall workflow is depicted in Figure 1A. Finally, we present the results in an accessible and searchable format and make it available on the web.

Adipocyte differentiation markers and analysis quality assessment
The mouse fibroblast cells, 3T3-L1 pre-adipocytes, differentiate into mature adipocytes when treated with the MDI differentiation medium (160 nM insulin, 250 nM dexamethasone, and 0.5 mM 1-methyl-3-isobutylxanthine). Once fully confluent, 3T3-L1 pre-adipocytes undergo a growth arrest followed by cell division and differentiation upon treatment with MDI [6]. Al Adhami et al. performed RNA sequencing of 8 samples of MDI-induced 3T3-L1 cells at 4 time points corresponding to these differentiation stages; day 0 (pre-adipocyte proliferation), day 2 (quiescence state), 10 h (clonal expansion) after MDI treatment and day 6 (differentiation) after MDI treatment [7]. To confirm the differentiation of the adipocytes in this dataset, we first examined the expression levels of 3 important adipocyte markers, which are correlated with the maturation of adipocytes as reported previously [8]; Cebpa, CCAAT/enhancer binding protein (C/EBP) alpha; Lpl, lipoprotein lipase; Pparg, peroxisome proliferator activated receptor gamma. The mRNA levels of these three adipocyte differentiation marker genes relatively increased in differentiating cells (day 6 after MDI treatment) compared to proliferating cells (day 0) ( Figure 1B). Additionally, several essential lipogenesis markers including Fasn, fatty acid synthase, and Acaca, acetyl-CoA carboxylase, were highly elevated in differentiating cells ( Figure 1C), indicating that adipogenesis-associated lipid biosynthesis is significantly activated.
To ensure the reliability of our analysis, we performed several steps to quality assess and validate the analysis. First, DESeq2 includes a per-gene estimate of dispersion, which is the within-group variability, in the modeling step. Supplementary Figure 1A shows the per-gene estimates, fitted and final estimates of dispersion versus the mean normalized counts for each gene. Second, the multi-dimensional scaling (MDS) analysis based on total (n = 19692) gene counts showed a clear separation between experimental groups (Supplementary Figure 1B). Finally, another dataset composed of 6 RNA-seq samples of MDI-induced 3T3-L1 at 2 time points (day 0 and day 7 after MDI treatment), deposited by Duteil et al. [9], was obtained and processed using the same pipeline. The average base expression of autophagy-related genes (n = 131) from both datasets showed a high correlation (Pearson's coefficient values = 0.87) between these two independent experiments (Supplementary Figure 1C). These suggest that the dataset we used here is suitable for the autophagy-dependent study in adipogenesis.

Gene expression of autophagy-related genes in adipogenesis
In order to examine how overall autophagy genes is regulated with respect to the differentiation stages of 3T3-L1 cells, we systematically analyzed the profiles of autophagy-related gene expression using the workflow model depicted in Figure 1A. First, we conducted a pairwise comparison of the autophagy-related genes expression among the 3 differentiation stages using the proliferating cells as a control, and the numerical summaries are shown in Figure 2A. Of total autophagyrelated genes (n = 131) identified in the gene ontology term (autophagy), 77 genes were expressed deferentially (adjusted p-value < 0.1) between quiescent (day 2) and proliferating cells (day 0), 93 genes between clonal expanding (10 h after MDI treatment) and proliferating cells (day 0), and finally 96 genes between differentiating (day 6 after MDI treatment) and proliferating cells (day 0). Using the Likelihood Ratio Test (LRT), we then tested the different expression pattern across the 4 differentiation stages. Several autophagy-related genes (n = 35) were found to be significant (adjusted p-value < 0.1 and log fold-change > 1) ( Figure 2C).

The distinct patterns of autophagy-related genes profiles during the adipocyte differentiation
In the Figure 2C, the heatmap suggests that there are 4 distinct patterns (temporal profiles) of autophagyrelated gene expression during the course of adipogenesis. To further characterize these patterns, we used the c-means algorithm to cluster the autophagy-related genes based on the Manhattan distance. Indeed, the resulting 4 clusters (cluster #1~4) corresponded to the main temporal profiles (Figure 3 top). Except the first cluster, autophagy-related genes exhibited relatively low expression at the proliferating stage (day 0) and a higher expression at the differentiating stage of adipocytes (day 6 after MDI treatment). While autophagy-related genes in cluster 3 showed a gradual increase of expression over time, clusters 2 and 4 showed that expression of some autophagy genes was occasionally decreased at 10 hours after MDI treatment. Genes in the cluster 1 showed a reverse pattern-higher expression at proliferating cells (day 0) and a lower expression at the differentiating cells (day 6 after MDI treatment). The members of each cluster are enumerated by symbol and gene name in Table 1.
Next, we used an algorithm to find other nonautophagy genes, not identified as autophagy-related genes in the gene ontology term (autophagy), with similar expression patterns. We first identified the clusters' mediods-a gene with the most typical pattern to a certain cluster-for each of the above described clusters and used the genefinder algorithm in the genefilter package to find the 5 closest genes to each of the clusters' mediods. We showed that the expression patterns of the closer 20 genes (4 groups: 5 genes in each group) to the previously identified clusters (Figure 3 bottom) and listed their symbols and gene names in Table 2.
In order to validate the expression pattern of genes belong to each autophagy cluster, we performed the differentiation of 3T3-L1 pre-adipocytes in the MDI differentiation medium as described in "Materials and Methods" and collected total RNAs from cells with four different stages; proliferation (day -2), quiescence (day 0), clonal expansion (10 h after MDI induction) and differentiation (6 days after MDI induction) shown in Figure 4A. The quantitative PCR data of some representative genes: Ubqln2 (cluster #1), Pink1  cluster #2), Foxo1 (cluster #3), and Prkaa2 (cluster #4), showed that the relative expression pattern of each gene at differentiation stages well agreed with the depicted patterns by members of the four autophagy clusters shown in Figure 3 and Table 1 ( Figure 4B).

Enrichment of autophagy gene sets during the adipocyte differentiation
To further understand the ongoing regulation of autophagy at the gene set level during adipogenesis, we first identified the autophagy offspring (sub-categories) gene sets in the gene ontology. Then, we used the gene count/ fragment per million (fpm) to test for their differential enrichment during courses of adipocyte differentiation. Few autophagy gene sets showed significant enrichment (False Discovery Rate < 0.1) between quiescent and proliferating cells. However, a high number of autophagy gene sets were significantly enriched between adipocyte differentiating cells and proliferating control cells due to a largely increased proportion of deferentially expressed genes between the two stages ( Figure 5). We present the proportion of deferentially expressed genes for each autophagy gene set at three different conditions (quiescence, clonal expansion, differentiation) Figure 2: Differential expression of autophagy-related genes and exons during adipocyte differentiation. Raw reads of 8 MDI-induced 3T3-L1 RNA-seq samples were aligned to the mm10 reference genome using HISAT2. Aligned reads were counted at feature; gene or exon level using HTSeq. Autophagy-related features were identified using the gene ontology term (autophagy) in the mouse organism annotation package org.Mm.eg.db. Feature counts were used to conduct pairwise comparisons between developmental stages (quiescence, clonal expansion and differentiation) and the control stage (proliferation). (A) A summary of DESeq2 differential expression analysis of (n = 131) autophagy-related genes. (B) A summary of DEXSeq differential exon usage of (n = 1777) autophagyrelated exons. Venn diagrams show significant features at (adjusted p-value < 0.1). (C) Log2 gene counts were used to conduct Likelihood Ration Test (LRT) for differential expression among developmental stages; proliferation, quiescence, clonal expansion and differentiation. 35 autophagy-related genes were found significant (adjusted p-value < 0.1 and log fold-change > 1) and their scaled counts are shown as a heatmap (dark green for high and light green for low counts). Rows are labeled by official symbols and columns by the SRA sample identifiers. Sample annotations (time point and induction) are mapped to the corresponding samples by a color key. www.oncotarget.com compared to the control condition (proliferation) in Figure 5 and listed the corresponding false discovery rates in Table 3. To better visualize the changes in the autophagy process during the adipocyte differentiation, we mapped the standardized average gene counts of two conditions (left, proliferation and right, differentiation) as color values to the mouse KEGG pathway for autophagy ( Figure 6).
Similarly, we identified the signaling pathways in the KEGG database that share at least one of their members with the autophagy regulation pathway. After testing for the enrichment of the pathways using the fragment per million gene counts, we listed the significant pathways in the later three differentiation stages (quiescence, clonal expansion, differentiation) compared to the earliest proliferation stage in Table 4. The mTOR, Jak-STAT, insulin and adipocytokine signaling pathways were highly enriched between differentiating and non-differentiating (proliferating) cells (Supplementary Figure 2), suggesting that the role of these signaling pathways in adipogenesis is associated, at least in part, with regulating autophagy.

Differential exon usage of autophagy genes in adipogenesis
In addition to the changes of the expression level of genes involved in the cellular process, relative usage of certain exons and the abundance of certain transcripts is an additional potential mechanism for regulation and modification of autophagy during the adipocyte differentiation. We tested for the differential exon usage in the eight samples of 3T3-L1 cells used in this study and identified potential alternative splicing events. A few of these events were identified in quiescent cells and more events at the clonal expansion and differentiating cells compared to the control proliferating cells. We summarized the differential exon usage events among autophagy-related exons (n = 1777) in Figure 2B. In addition, we presented the relative usage of 15 exons of Acbd5 (acyl-CoA binding domain containing 5), a gene known to function in the transport and distribution of long chain acyl-Coenzyme A in cells, during the adipocyte differentiation in Figure 7. Exons 1, 2 and 8

Figure 3: Expression patterns of autophagy-related genes and non-autophagy genes with similar patterns. (top)
Autophagy-related genes were identified using the gene ontology term (autophagy) in the mouse organism annotation package org.Mm.eg. db., then were clustered using the Manhattan distance and the c-means algorithm. Each panel represents a distinct cluster (expression profile over time); gray color indicates the area spanned by the cluster members, and the blue line is the linear trend within the cluster.
y-axis represents the standardized expression (mean subtracted and divided by the standard deviation). (bottom) Non-autophagy genes were examined for showing expression patterns similar to the 4 autophagy-related genes' clusters using the genefinder algorithm in the genefilter package. The 5 closest genes (indicated by different color lines) within each cluster were selected and shown in separate panels. y-axis represents the standardized expression (mean subtracted and divided by the standard deviation). showed significant relative usage over time (adjusted p-value < 0.1). Exon 2 is only a part of 2 transcripts (bottom), which indicates a potential splicing event. In other words, transcripts containing exon 2 were less expressed in differentiating cells compared to earlier stages and probably the first two transcripts were more abundant.

DISCUSSION
Autophagy is an essential process during the preadipocyte 3T3-L1 differentiation into mature adipocytes. According to previous studies, inhibition of autophagy by knockdown of Atg5 or Atg7 genes in 3T3-L1 pre-adipocytes prevents the adipogenesis and the formation of lipid droplets. Also, autophagy inhibition causes a significant decrease in adipocyte differentiation markers, such as CCAAT/enhancer-binding protein-alpha (CEBPa) and peroxisome proliferator activated receptor-gamma (PPARg) [10]. In this report, we examined the differential expression of autophagy-related genes at the transcriptional level during the adipocyte differentiation using public-access RNA-seq datasets. Our analysis reveals that the expression of Atg5 and Atg7 genes is maintained at a high level during the course of differentiation with respect to the maturation   stage when 3T3-L1 pre-adipocytes are induced in the differentiation medium (data not shown). Although we cannot quantify the lipid contents in this study, we show the relatively high expression of a number of protein markers that are highly correlated with lipogenesis in adipose tissue in Figure 1C [11].
Other cells than adipocytes also have the ability to store different amounts of lipids, for instance, hepatocytes and neurons; in some pathological conditions [10,12]. In particular, hepatocytes use autophagy to generate free lipids from lipid droplets for energy metabolism in a specialized lipid degradation process called lipophagy [10]. In addition to the role of autophagy in lipid droplet formation in adipogenesis, autophagy also takes a significant part in the lipolytic pathway. The gene ontology term lipophagy represented by one gene member, RAS oncogene family (Rab7), appears to be slightly activated in differentiating adipocyte compared to the control pre-adipocytes. In addition, adrenergic receptor, beta 2 (Adrb2), a member of the term positive regulation of lipophagy is significantly up-regulated in the same comparison. However, how these two lipophagy-associated genes at the gene ontology term can function in the canonical autophagy process still remains to be demonstrated.
The classic ability of autophagy to recycle the cellular components is particularly important during the adipocyte differentiation. Autophagy targets different intracellular organelles, cytoplasmic regions and protein aggregates. This in turn allows adipocytes to store large amounts of lipids and to gain a distinct morphology. Additionally, the nutrients acquired by the autophagy-dependent degradation in adipogenesis can be good sources for adipocyte growth and lipid biosynthesis. The gene set enrichment analysis reveals activation of autophagy gene subsets responding to distinct intracellular targets; mitophagy (mitochondria), nucleophagy (nucleus), reticulophgay (endoplasmic reticulum), pexophagy (peroxisomes), xenophagy (cytoplasmic regions) and aggrephagy (protein aggregates). Moreover, mitophagy-specific genes such as Pink1, or Park2 [13][14][15] are differently expressed at different time Figure 5: Enrichment of autophagy gene sets during adipocyte differentiation. Autophagy-related genes were identified using the gene ontology term (autophagy) in the mouse organism annotation package org.Mm.eg.db. Package limma were used to test for the gene set enrichment by applying the rotation test to the fragment per million (fpm) matrix and the comparisons between proliferating cells (day 0) and quiescent cells (day 2), clonal expanding cells (10 h after MDI treatment) and differentiating cells (day 6 after MDI treatment). Bars represent the proportion of genes in a gene set (labeled by GO ids) that are deferentially expressed at a certain comparison. Color indicates whether the overall gene set enrichment is significant (false discovery rate < 0.1) (green, yes or red, not). www.oncotarget.com  [16,17].
Optn gene also plays a critical role in degradation of some cytoplasmic regions via xenophagy [18]. Indeed, we showed that the expression of these two genes was regulated in a developmental stage-specific manner during adipocyte differentiation. Given the role of autophagy in adipocyte differentiation for lipid droplet formation and recycling the cellular components, we might also expect some significant quantitative changes in the process or in its associated regulatory pathways. It is not clear whether the quantification of autophagosomes at a certain time point without blocking the autophagy flux is possible. Nonetheless, transcript levels of both microtubule associated protein 1 light chain 3 (maplc3a) and autophagy-related 5 (Atg5) are highly correlated to autophagosome formation. As shown in ( Figure 5 and Table 3), we observed a significant up-regulation of members of the autophagosome assembly genes (31% and 75% up-regulation at 10 h and 6 days after induction, respectively; adjusted p-value < .01) Consequently, we also observe a significant enrichment of autophagy-related gene ontology terms; namely protein localization to pre-autophagosomal structure and protein lipidation involved in autophagosome assembly. Taken together, we can infer a quantitative increase of autophagy during the adipocyte differentiation based on increased transcription levels of Map1lc3a and other genes involved in its localization and lipidation.
Sub-grouping process of autophagy to several offspring gene sets in the gene ontology can give the detailed information to understand the role of autophagy in adipogenesis at the gene set level. Based on only mRNA levels obtained from public-access RNA-seq data, we indeed found a significant enrichment of many autophagy gene sets during the course of adipocyte differentiation although there still would be difficult to understand how each autophagy gene set is exactly related to among different types of autophagy.
Detailed knowledge of relative gene expression levels of key genes in autophagy-related gene sets and their overall enrichment over time could be used to emphasis certain aspects of the involvement of autophagy in adipocyte differentiation. Basal autophagy is activated when pre-adipocytes respond to the differentiation stimulus. Organelle-specific autophagy is probably employed to achieve the structural modification and meet the metabolic demands in mature adipocytes. Quantitative and regulatory changes are indeed required to govern the autophagy response during differentiation. Not only autophagy-related gene sets such as autophagosome assembly and signaling pathways such as mTOR and Jak-STAT but also adipogenic pathways like insulin and adipocytokine signaling show a significant enrichment at the differentiation stage of adipocytes compared to the control stage.
Much of the regulation of autophagy in adipogenesis is probably happening at protein and protein modification levels, and wouldn't be easily captured by the changes at the mRNA level. Using parallel datasets from proteomic, RNA-seq or transcription factors ChIP-seq studies is one way to overcome the limitation of general protein functional approaches shown in many other studies [2,3,10]. Also, assessing the lipid contents or the autophagosome formation in a certain condition is not feasible by sequencing data alone. Therefore, integrating datasets of varies measurement technologies is a plausible future application of these findings.
In conclusion, the systematic analysis using public-access databases such as microarray or RNA-seq can provide new valuable information to collectively understand a specific cellular pathway. Here, we show that expression of many autophagy-related genes are highly up-regulated at the differentiation stages of adipocytes in accordance with previous reports, and autophagy gene sets defined in the gene ontology are enriched among the developmental stages, further confirming the critical role of autophagy in adipogenesis.

Raw data and processing
The dataset (Al Adhami et al.) consists of 8 samples (two of each time point) of 3T3-L1 pre-adipocytes treated with/without MDI (160 nM insulin, 250 nM Figure 6: Autophagy KEGG pathway map in proliferating pre-adipocytes and differentiating adipocytes. Mouse autophagy pathway (04140) graph was obtained from the KEGG pathway database using pathview package. Scaled fragment per million (fpm) of genes were used to map the expression levels (red for high and green for low) to the graph. Each node is divided to two parts; left, proliferating pre-adipocytes (day 0) and right, differentiating adipocyte cells (day 6 after MDI treatment). .01 False discovery rates of rotation tests of KEGG pathways that shares one or more of autophagy-related genes in three differentiation stages compared to control dexamethasone, and 0.5 mM 1-methyl-3-isobutylxanthine) at four time points; day 0 (pre-adipocyte proliferation), day 2 (quiescence state) and 10 h (clonal expansion) after MDI treatment, day 6 (differentiation) after MDI treatment. The protocol for generating the RNA libraries and sequencing is described at [7]. Raw reads for 8 samples were obtained from ENA (PRJNA218101). Reads were aligned to the mouse UCSC (mm10) reference genome using the splicing aligner HISAT2 [19]. Aligned reads were counted using HTSeq [20] and summarized at gene and exon level guided by the corresponding annotation data from the UCSC mm10. A validation dataset (Duteil et al) consists of 6 samples (three of each time point) of 3T3-L1 cells at two time points day 0 and day 7 after MDI treatment. The protocol for generating the RNA libraries and sequencing is described at [9]. Raw data were obtained from ENA (PRJNA219405) and processed by the same above mentioned pipeline to summarize counts at gene level.

Differential expression and exon usage
Pairwise comparison between control (proliferation) and 3 conditions (quiescence, clonal expansion and differentiation) were conducted by applying Wald test of the negative binomial distribution to the log2 gene counts using Figure 7: Differential exon usage of Acbd5. Counts of total 16073 exons in 8 samples were used to test for differential usage using DEXSeq for the pairwise comparison between proliferating samples (day 0) and quiescence (day 2), clonal expansion (10 h after MDI treatment) and differentiating samples (day 6 after MDI treatment). Exon usage values of 19 exons of Acbd5 at different differentiation stages (green, proliferation; purple, quiescence; red, clonal expansion; blue, differentiation) are shown with their corresponding transcripts and genomic ranges at the bottom. www.oncotarget.com DESeq2 [21]. Likelihood Ratio Test (LRT) was also applied using DESeq2 to test for differential expression across all conditions. DEXSeq was used to test for the differential exon usage [22]. p-values were adjusted for multiple testing using BH and the cutoff 0.1 was chosen. DESeq2 and DEXSeq estimate a size factor, that is the median of the ratios of the observed counts, one value for each sample as a normalization factor to correct for potential sequencing depth and other biases to make gene measurements form different samples comparable. We limited the results to the group of genes of known function or relation to autophagy as identified by the gene ontology term autophagy (GO:0006914) and referred to them as autophagy-related genes.

Clustering and pattern finding
C-means algorithm in the e1071 package was used to cluster the genes into distinct temporal profile. C-mean is a fuzzy version of the k-means algorithm [23]. Package genefilter was used to find genes with similar expression pattern to the clusters' mediods-genes with the most typical pattern to identified cluster [24]. For a certain cluster, genefinder calculates the Euclidean distance of all genes to the a mediod, ranks genes by their distance and chooses the n closest genes to the mediod.

Gene set enrichment analysis
Gene ontology (GO) was used to identify the autophagy gene sets and related genes [25]. KEGG pathways were used to identify signaling pathways that share one or more of the autophagy related genes [26]. The mouse organism annotation package (org.Mm.eg.db) was used to retrieve the relative annotations [27]. Gene counts were first transformed to the fragment per million and using package limma a rotation test for the gene set and pathway enrichment was applied to the comparison between control (proliferation, day 0) and 3 conditions (quiescence, day 2; clonal expansion, 10 h after treatment; differentiation, day 6 after treatment) [28]. p-values were adjusted for multiple testing using False Discovery Rate (FDR) and the cutoff 0.1.

Cell culture, MDI induction, and qPCR
3T3-L1 pre-adipocytes were cultured and maintained in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% calf serum (Gibco, Carlsbad, CA, USA) at 37° C in a humidified atmosphere of 5% CO 2 . At 100% confluence, cells were transferred to a fresh medium and incubated for 2 days. Arrested cell were then moved to fresh DMEM supplemented with 10% FBS and the induction mixture MDI (160 nM insulin, 250 nM dexamethasone and 0.5 nM 1-methyl-3isobutylxanthine). The medium was then changed twice at 2 days' interval to DMEM + 10% FBS and in addition to 160 nM at the first change only. Cells were harvested at 4 time points corresponding to the 4 differentiation stages; proliferation (day -2), quiescence (day 0), clonal expansion (10 hours after MDI induction), and differentiation (6 days after MDI induction). Total RNAs were purified using TRIzol and assessed by measuring the 260/230 diffraction ratios to ensure the quality of RNA samples and to calculate the proper amount for the qPCR quantification.
One-step real-time quantitative PCR (QuantiNova SYBR Green PCR kit, QIAGEN, Rotor-Gene Q) was used to amplify and quantify the mRNA levels in three independent samples from the 4 time points. The primers used for qPCR were as follows; mUbqln2 (forward-TTGAGCTGTTCCAGTTGCTG, reverse-ACCCAA CCAGCAGTTCATTC), mPink1 (forward-AGTGTCC AGTGGGTCAGACA, reverse-CTGATCGAGGAGAA GCAGG), mFoxo1 (forward-TGCTGTGAAGGGACA GATTG, reverse-GAGTGGATGGTGAAGAGCGT), and mPrkaa2 (forward-ACTGCCACTTTATGGCCTGT, reverse-GATCGGACACTACGTCCTGG), and m18S (forward-ACCGCAGCTAGGAATAATGGA, reverse-GCCTCAGTTCCGAAAACCA). The pcr package was used to apply the comparative CT method on the data. CT values were normalized to that of the control 18S and calibrated to the first time point (proliferation stage).

Software environment and reproducibility
The data analysis and visualization were mainly done in an R environment. A fully reproducible version of this manuscript is available at (https://github.com/ MahShaaban/3T3-L1) along with the scripts that were used to process and analyze the data.