Mouse Embryonic Fibroblast Adipogenic Potential: A Comprehensive Transcriptome Analysis

ABSTRACT Our understanding of adipose tissue has progressed from an inert tissue for energy storage to be one of the largest endocrine organs regulating metabolic homoeostasis through its ability to synthesize and release various adipokines that regulate a myriad of pathways. The field of adipose tissue biology is growing due to this association with various chronic metabolic diseases. An important process in the regulation of adipose tissue biology is adipogenesis, which is the formation of new adipocytes. Investigating adipogenesis in vitro is currently a focus for identifying factors that might be utilized in clinically. A powerful tool for such work is high-throughput sequencing which can rapidly identify changes at gene expression level. Various cell models exist for studying adipogenesis and has been used in high-throughput studies, yet little is known about transcriptome profile that underlies adipogenesis in mouse embryonic fibroblasts. This study utilizes RNA-sequencing and computational analysis with DESeq2, gene ontology, protein–protein networks, and robust rank analysis to understand adipogenesis in mouse embryonic fibroblasts in-depth. Our analyses confirmed the requirement of mitotic clonal expansion prior to adipogenesis in this cell model and highlight the role of Cebpa and Cebpb in regulating adipogenesis through interactions of large numbers of genes.


Introduction
Adipose Tissue (AT) function has advanced from its role as a primary storage of triglyceride in the form fat droplets for energy to the secretion of various signalling factors, known as adipokines that regulates various metabolic pathways. Thus, AT has been characterized as a major endocrine organ [1]. Adipogenesis, the transformation of preadipocytes into mature adipocytes, plays a critical role in development of specific functions and differentiation of cellular subtypes [2][3][4].
Adipocyte differentiation is a dynamic process tightly regulated through a multistep process conserved across species [2]. The process involves temporal expression of various signalling cascades that regulate the expression of transcription factors and thereby proadipogenic genes such as peroxisome proliferatoractivated receptor-у (PPARG) [5,6]. Adipocyte differentiation is characterized by significant changes to cell morphology -from a fibroblast to a rounded spherical form -and acquisition of functional characteristics of AT cells [7]. These properties are directed by the expression of PPARG, which is central to adipogenesis. This receptor is encoded by Pparg and cooperate with members of the CCAAT/enhancer-binding proteins (CEBPs) family of transcription factors, specifically CEBP-β (CEBPB), encoded by Cebpb, and CEBP-α (CEBPA), encoded by Cebpa. Several studies have demonstrated the interplay among these factors during adipocyte differentiation and in enhancing adipogenic cell fate [8][9][10][11]. Mechanistically, induction of CEBPB, at an early phase of adipogenesis, promotes induction of PPARG and CEBPA by binding to their proximal promoter regions. Therefore, this results in activating several pro-adipogenic genes [12][13][14][15][16] that support the progression of preadipocytes cells into mature adipocytes [11,17].
Global gene profile changes in various murine cell models are reported in the literature that show the biological and clinical significance of adipogenesis through comprehensive analysis of differentially expressed genes (DEGs) using tools such as RNA sequencing (RNA-Seq) [18][19][20][21][22][23][24].
Through statistical significance and fold induction, the final lists of DEGs are further analysed using gene ontology (GO) analysis, and can assist in construction of proteinprotein interactions (PPI) and gene networks. However, RNA-Seq data are variable, and biased analyses are possible. Thus, unbiased methods are necessary for generating transcriptome data to increase confidence of DEG lists. The robust rank aggregation (RRA) method that relies on algorithms that produce a significant score for genes that consistently rank better than expected has been applied to transcriptome data [25,26]. In this method estimates of significance are provided as p-values that assign probabilities showing that the differential expression of a given gene is related to adipogenesis [25,27]. p-values can be used to rank DEGs, where the lower the value, the higher its rank. An RRA approach allows multiple transcriptome profiles to be characterized by more reliable molecular targets.
In the current study, adipogenesis in MEFs transcriptome profiles was explored using generated DEGs at different timepoints. DEG lists relevant to adipose were generated based on statistical significance and used for an RRA method, and in construction of PPI and gene networks. Therefore, this study aims to utilizes transcriptome analysis methodology to acquire resourceful understanding in gene expression changes in MEFs that underlie adipocyte differentiation.

Oil red O staining and quantification
Differentiated and undifferentiated MEFs were washed with 1 ml of 1X PBS (Corning) once before fixing for 5 min with 4% paraformaldehyde in PBS at room temperature. After fixation, MEFs were washed three times with 1 ml of PBS and once with 60% isopropanol and were completely air-dried. MEFs were then stained with Oil Red O (ORO) [6 ORO:4 ddH 2 O ratio] (Sigma©) for 10 min at room temperature. Excess stain residue was removed with four ddH 2 O washes. Approximately 1 ml of PBS was then added for microscopic visualization. Images were processed at 460x magnification using a 20x objective lens, with transmitted bright-field light (EVOS® FLoid® Cell Imaging Station). For quantification, ORO stain particles were eluted with 100% isopropanol and analysed using Thermo Scientific™ Varioskan® Flash for spectrophotometry readings at 514 nm.

RNA isolation and purification
Total RNA was extracted from MEFs using a combination of TriZol (Life Technologies™) and RNeasy Mini Kit (Qiagen©), with modifications from the manufacturer's protocols. Cells in each well were washed once with 1 ml of PBS prior to the addition of 1 ml of TriZol reagent. TriZol lysates were added to fresh 1.5-ml tubes. Approximately 0.2 ml of chloroform was added per ml of TriZol and was centrifuged at a speed of 12,000 rpm at 4°C for 15 min. The upper aqueous phase (~400 ml in volume) was transferred to a fresh 1.5-ml tube. An equivalent volume of 70% ethanol was mixed with the RNA solution, and downstream purification was performed using a RNeasy Mini Spin Column (Qiagen©) according to the manufacturer's protocol. RNA concentrations were determined using a Thermo Scientific™ NanoDrop 2000™.

cDNA synthesis and quantitative PCR
Total isolated RNA (0.1 ng) was used as a template for synthesizing complementary DNA (cDNA) using a First-Strand Synthesis Kit (Invitrogen™) according to the manufacturer's protocol. cDNA was diluted 1:50 in ddH 2 O prior to its use as a template for real-time PCR (RT-PCR) using SYBR Green (Thermo Scientific™). A total volume of 10 ul was prepared for the SYBR Green qPCR reaction [2.5 ul diluted cDNA product or nH 2 O, 2.5 ul 2 mM forward + reverse primer mix (Integrated DNA Technologies, Inc.) (Table 8) and 5 ul of SYBR Green]. The qPCR thermal cycle utilized an Applied Biosystems™ StepOnePlus™ platform, with one cycle of 50°C for 2 min and 95°C for 10 min followed by 44 cycles of 95°C for 15 s and 60°C for 1 min. Expression of gene targets values was normalized to housekeeping gene β-actin and was estimated using the ΔΔCt approach. Results represented on a log2 scale.

RNA-Seq library preparation and sequencing
Total RNA quality was estimated based on RNA integrity number (RIN) using Agilent BioAnalyzer 2100™. RNA samples with a RIN>8 was used for library preparation. RNA-Seq libraries were prepared using an Illumina® TruSeq Stranded mRNA Prep Kit accordingly with the manufacturer's LS protocol. Samples were barcoded, multiplexed and sequenced (100 bp pair end) using the Illumina® NextSeq 550 platform at NYU Abu Dhabi (NYUAD) Genomic Core facility (Abu Dhabi, U.A.E).

Transcriptome data generation
The DESeq2 computational pipeline was used to estimate raw count reads aligned to the reference genome [35]. Mouse genome (GRCm38/mm10) from the University of California Santa Cruz Genome Browser (https://genome.ucsc.edu/) was utilized as a reference [36]. Computing methods used a Linux-based command system on the New York University Abu Dhabi (NYUAD) high-performance computing (HPC) server platform Dalma (https://wikis.nyu.edu/display/ADRC/ Cluster+-+Dalma). Raw data counts were deposited at National Centre for Biotechnology Information (NCBI) database and are available under Gene Expression Omnibus (GEO) accession number GSE152750.

Bioinformatic and computational analysis
Heatmaps and correlation through PCA and distance heatmap dendrogram analysis were generated by RNA-Seq START (Shiny Transcriptome Analysis Resource Tool), via the New York University Abu Dhabi Centre of Genomic and Systems Biology (NYUAD-CGSB) Bioinformatics Online Analysis and Visualization Portal (http://tsar.abudhabi.nyu.edu/) [37].
A PPI network of differentially expressed genes (DEGs) was constructed with the STRING database (http://www. string-db.org) followed by visualization with the Cytoscape software (version 3.6.0, Washington, DC, USA) (http:// www.cytoscape.org/). The degree of nodes was determined with the plug-in Network Analysis app on the Cytoscape software. Gene interaction network models were constructed with the Biological General Repository for Interaction Datasets (BIOGRID; https://www.thebiogrid. org/) database followed by visualization of DEGs with the Cytoscape software.

Statistical analysis
Statistical analysis used Student's t-tests in Microsoft Excel™ (Microsoft®). Results are represented as the mean of at least three independent experiments, and a p-value<0.05 was considered significant.

Formation of adipocyte-like cells from MEFs
We used MEFs in an initial incubation of cellsfrom day 0 (D0) -with adipogenic differentiating factors in culture medium (see Material and Method section) for 3 days (D3) followed by the removal of these factors and an additional incubation for 2 days (D5) (Figure 1a). Lipid accumulation was assessed with Oil Red O (ORO) staining on days 3 and 5. As expected, lipid accumulation was observed as droplet formation in treated (+) cells, but not in untreated (-) cells at both timepoints. Qualitatively, more and larger droplets were seen on D5 (Figure 1b). This finding was confirmed quantitatively using spectrophotometry with extracted stain. Higher optical density (OD) was recorded incrementally from treated cells, where OD was higher on D5 than on D3 (Figure 1c). We also examined endogenous expression of specific gene markers that are constitutively active in AT to confirm that results were not due to a random effect of treatment. Expression of mRNA levels for the adipocyte protein 2 (aP2) gene, also known as free fatty acid-binding protein 4 (Fapb4) [40] (Figure 1d; left chart), lipoprotein lipase (Lpl) gene [41] (Figure 1d; middle chart), and adiponectin (Adipoq) gene [42] (Figure 1d; right chart) were significantly upregulated after treatment. Comparatively, mRNA expression levels of the aP2, Lpl, and Adipoq genes were higher on D3, but not significantly, than on D5. These observations indicate that MEF developed adipogenic identity and functional activity. Thus, the process used to induce cells towards differentiation into adipocytes is a useful platform for transcriptome analysis.
Based on the correlation analysis, DEGs were identified for intermediate (MEFD3plu/MEFD3min) and terminal timepoints (MEFD5plu/MEFD5min). DEGs were genes that showed a statistically significant change in expression (p-adj<0.05) between treatments. Using volcano plots, a total of 9186 ( Figure 2e) and 8226 (figure 2f) DEGs were identified for MEFD3plu/ MEFD3min and MEFD5plu/MEFD5min, respectively. Differences in numbers of DEGs may reflect a series of synergetic processes, such as genes response to signalling molecules (i.e. cellular dynamics) or nuclear changes (i.e. methylation/acetylation events) during adipogenesis. Such processes may be involved in specifying cell fate or enhancing adipogenic features. DEGs from both timepoints were further analysed via log scaled variance-stabilizing transformation (VST) to generate a heatmap of hierarchal clustering. Expression of selected DEGs lists in each timepoint was analysed across all conditions (MEFD0, MEFD3min, MEFD3plu, MEFD5min, and MEFD5plu) and their corresponding replicates. Variable expression patterns were observed across conditions (Figure 2c and D), which may reflect potency of the differentiation treatment. However, between both timepoints, there were no significant changes in the expression patterns and number of DEGs list isolated from MEFD3plu/MEFD3min and MEFD5plu/MEFD5min.

Gene ontology, expression, and network interaction analysis of genes related to fat differentiation and development
Functional aspects of DEGs were examined using gene GO databases. Firstly, DEGs from both timepoint sets were organized into a Venn diagram to eliminate redundancy of genes ( Figure 3a). The diagram showed 5314 DEGs shared between MEFD3plu/MEFD3min and MEFD5plu/MEFD5min and 3872 and 2912 uniquely expressed genes between MEFD3plu/ MEFD3min and MEFD5plu/MEFD5min, respectively.
The common 5314 DEGs showed 3224 enrichment terms overall. The enrichment list showed 92 enrichment terms related to differentiation processes, including neuron, chondrocyte, osteoblast, macrophage, myeloid, and various immune and progenitor cells. Interestingly, among enrichments, GO terms fat cell differentiation (GO:0045444, 108 genes) and brown fat cell differentiation (GO:0050873, 32 genes) were most significantly present. Further filtering of enrichment terms was performed with focused terms relevant to fat cells and adipose tissue ( Table 1) Table 1). MEFD3plu/MEF3min unique DEGs also relevant enrichments towards fat development including regulation of fat cell differentiation (GO:0045598,21 genes; Table S1), positive regulation of fat cell differentiation (GO:0045600, 11 genes; Table S2), and fat cell differentiation (GO:0045444, 26 genes; Table S3) (lower panel: Figure 3b). Whereas MEFD5plu/MEFD5min unique DEGs, constrictively showed enrichments towards brown fat cells including positive regulation of brown fat cell differentiation (GO:0090336, 2 genes; Table S4), regulation of brown fat cell differentiation (GO:0090335, 2 genes; Table S5), and brown fat cell differentiation (GO:0050873, 2 genes; Table S6) (middle panel: Figure 3b). Given that the common DEGs showed relevant enrichment terms as opposed to unique DEGs. Thus, the former was further analysed.
Expression patterns of the genes identified in the common enrichment list were further analysed using a metric heatmap (Supplementary Figure 2). The heatmap included all DEGs in each previous enrichment into a single list (119 genes), by removing redundant genes and normalizing gene counts on a log2 scale. The identified genes were categorized according to the timing of up or down differential expression (Supplementary Figure  2  Specific interactions among these genes may be necessary for any given biological process. An interaction network is essential for understanding time-dependent regulation during differentiation observed in the previous heatmap with some genes. Thus, the relationship between gene expression patterns and timing is essential for elucidating fat differentiation and adipose development. Furthermore, key genes may be primary controllers of regulatory pathways. To explore this, a PPI network was constructed and curated via STRING and Cystoscape [46,47], respectively (Supplementary Figure 3). The numbers of degrees (adjacent links) from each node (gene) intersection were systematically selected with betweenness and closeness centrality ( Figure 4a, Table 2 -displaying the top 30 genes). Gene Akt1 showed the highest number of degrees, at 39, among all genes. Adipoq and Fabp4 (aP2) showed 21 and 16 degrees, respectively, and Ucp1 and Lpl showed 17 and 14 degrees, respectively. The top 30 genes from the GO term were further analysed at the biological process level with GO term analysis (Figure 4b) where most were associated with fat cell differentiation and metabolic and organic response processes (Table 3).

Gene interaction profile changes related to Cebpa and Cebpb expression
Most genes identified in the PPI network are known to influence adipogenesis, and Cebpa and Cebpb were further investigated for their relation to protein factors. The endogenous expression of these factors may differ in terms of interaction with other proteins in a timedependent manner (MEFD3plu/MEFD3min and MEFD5plu/MEFD5min). Therefore, Cebpa and Cebpb were examined with RT-PCR to determine and validate expression at the mRNA level (Figure 5a). On D3, significant upregulation of Cebpa and Cebpb was observed; however, on D5, only Cebpa showed significant upregulation (Figure 5a). These observations reflect the sensitivity of RNA-Seq analysis and suggest time-dependent interactions during treatment.
Both Cebpa and Cebpb are noted for their pioneering effects in driving adipogenesis. Therefore, gene interaction networks were formulated for both genes and were comparatively analysed against DEGs list from each MEFD3plu/MEFD3min and MEFD5plu/MEFD5min lists. Gene interaction networks for murine sources were downloaded from BioGRID and uploaded to Cytoscape [47,48]. The Cebpa network, consisting of 117 interactions, 63 and 35 genes were observed in DEGs list of MEFD3plu/MEFD3min (upper panel: Figure 5b; Table 4) and MEFD5plu/MEFD5min (upper panel: Figure 5c; Table 6), respectively. Whereas in the Cebpb network, consisting of 16 interactions, 11 and 8 genes were observed in DEGs list of MEFD3plu/MEFD3min (lower panel: Figure 5b; Table  5) and MEFD5plu/MEFD5min (lower panel: Figure 4c; Table 7), respectively. Most interactions identified are directed through a protein targeted manner.
To further decipher Cebpa and Cebpb interactions, PPI networks of the DEGs detected in each timepoint were constructed using STRING. Design of construct, formulated via Cytoscape, were restricted on direct interactions of DEG nodes with either Cebpa and Cebpb in each timepoint. The genes were ranked in degree of interaction in each timepoint after merging Cebpa and Cebpb profiles (left panels: Figure 6a and B). Based on the network maps, there were more direct interactions at day 3 (8 ranked orders) than that of day 5 (4 ranked orders). At day 3, Cebpa showed the highest level of interaction at 12 degree; whereas Cebpb was ranked in the fourth order, along with Cdk1 and Psmb10, at 5 degrees (right panel: Figure 6a). Interestingly, Hdac1 and Cdk4 were ranked in the second and third order at 9 and 6 degrees, respectively. Moreover, Wrd5 was ranked at the fifth order at 4 degrees; whereas Smarca2, Trib1, Rb1, Mapk14, and Gsk3b were collectively ranked in the sixth order with 3 degrees each. Further to those, Jpr2, Runx1t1, and Med1 were also collectively ranked at seventh order at 2 degree; whereas Pteg2 and Rnf1 were ranked at the eighth order with a degree level of 1. Inversely, at day 3, Cebpb showed the highest level of interaction at 6 degrees followed by Cebpa at 4 degrees (right panel: Figure 6b). Genes Bhlhe41 and Gsk3b were ranked at third order at 3 degrees each; whereas Cdk2, Sin3a, and Kmt2d were ranked at the fourth order with 2 degrees each. The expression of the genes was further analysed using a metric heatmaps, which showed variable patterns among the genes at different timepoints (MEFD0 -day 0, MEFD3minuntreated day 3, MEFD5min -untreated day 5, MEFD3plu -treated day 3, and MEFD5plu -  (Figure 6c). As previously shown, Cebpa, Cebpb, and Ptge2 were observed to be differentially upregulated in the treated conditions of day 3 and 5.
Interesting, Psmb10 and Hdac1 were observed to be differentially upregulated in treated MEF at day 5 as opposed to that at day 3. Whereas, Mapk14, Gsk3b, Trib1, Kmt2d, and Rb1 were differentially upregulated in treated MEFs at day 3 rather than day 5. There were also number of genes that were observed to be differentially downregulated due to the potency effect of differentiation treated, which included Cdk2, Med19, Cdk4, Jdp2, Rnf41, Ncoa3, Wdr5, Sin3a, Cdk1, and Bhlhe41. Overall, these variable observations between the different timepoints may suggest that specific interactions occur in a timedependent manner for Cebpa and Cebpb.

Discussion
AT consists of various cell types, with adipocytes as the major cell type present in mature tissues. Progenitor populations of stem cells and preadipocytes are also present. These latter cells can differentiate into various cell types, including adipocytes and its corresponding subtypes. Proliferation and differentiation of preadipocyte to mature cells continually sustain and maintain AT function. AT is present throughout the life cycle, and 50% of human subcutaneous fat cells are renewed every 8 years [49], illustrating the necessity of dynamic adipocyte replenishment through adipogenesis. Furthermore, from a global public health perspective, dysfunction of AT is linked to chronic diseases/disorder (i.e. obesity, cardiovascular diseases, diabetes) [4]. Therefore, the study of adipogenesis in vitro at molecular and mechanistic levels is essential for development of better therapeutic targets. Many studies that elucidate the processes of adipogenesis have used different cell lines, such as 3T3-L1, 3T3-F442A, MSCs, murine and human adipose stem cells (ASCs), and MEFs [16,[28][29][30][31][32][33][34]. Some of these same models have also been used to explore transcriptomic aspects of adipogenesis [16,19,20,[50][51][52]. Though MEFs have been used extensively, yet little is known about the dynamics of adipogenesis at a global gene level. MEFs are fibroblasts derived and isolated from mouse embryos. When cultured in vitro, MEFs display spindlelike features that are typical of many fibroblasts [53][54][55].
MEFs are used extensively as a cellular 'scaffolding' tool in stem cell biology in maintaining the growth of primitive undifferentiated cells such as mouse-and human-induced pluripotent or embryonic stem cells [55,56]. Moreover, in addition to adipocytes, MEFs display multilineage differentiation ability for the production of effector cell types, including osteocytes, chondrocytes and neurons [32-34,-57-59]. Most studies of adipocytes illustrate MEF adipogenic potential through observation of cellular morphological changes, such as lipid accumulation during differentiation, or by evaluating endogenous expression of known gene markers. In this study, MEFs were also observed to develop morphological changes associated with adipogenesis. MEFs accumulated lipids, shown by ORO staining, known to target lipid triglyceride droplets [28,29,60]. This accumulation was measured quantitatively with spectrophotometry [60]. Temporal lipid accumulation was assessed at 0, 3, and 5 days in cells in non-treated, basal medium, and differentiation medium groups. These time points were chosen due to sufficiency of lipid formation in MEFs when exposed to   16 26 differentiation factors, as observed previously [16]. Where continues incubation of MEFs with differentiation factors can lead continuous cell death [16], which will have biases towards exploring RNA-Seq data especially when targeting DEGs related to adipocyte or fat formation. Thus, a beginning (day 0), intermediate (day 3), and a terminal (day 5) time point was designated for this study. Further to this, the change from differentiation media to basal media was used to revive existing formed fat cells and explored its genotypic nature via transcriptome analysis. The differentiation media used in this study contained known inducing factors, including dexamethasone, insulin and 3-isobutyl-1-methylxanthine (Figure 1a). Although non-homogenous in lipid content, a positive linear increment in lipid droplets was seen in treated as opposed to non-treated cells. In parallel, similar increments were observed at the gene level with the expression of a lipolytic protein secreted during adipogenesis that is an in vivo feature of adipocytes including aP2, Lpl, and Adipoq [41,42,61]. Consistent with other Table 4. List of DEGs observed gene interaction network of Cebpa at MEFD3plu/MEFD3min highlighting interaction type, expression fold change (*) and its significance (**), and significance of interaction (***). studies, the present study observed formation of adipocyte-like features in vitro from MEFs under exposure of enhancing adipogenic factors. However, due to heterogenous population observed in MEFs, these observations do suggest that 3T3-L1 still remains to be an optimum model to study the functional aspect of Table 5. List of DEGs observed gene interaction network of Cebpb at MEFD3plu/MEFD3min highlighting interaction type, expression fold change (*) and its significance (**), and significance of interaction (***).  Table 6. List of DEGs observed gene interaction network of Cebpa at MEFD5plu/MEFD5min highlighting interaction type, expression fold change (*) and its significance (**), and significance of interaction (***).  adipogenesis and its metabolism. Yet, MEFs use as a model in AT development can be further explored and potentially utilized to study adipogenesis.
The ultimate aim of this study was to explore global gene expression controlled in vitro experiments. The study included analysis of transcriptome data based on statistical significance and ranking methods (RRA). Analysis of collected data offers valuable insight, biologically and clinically, to adipogenesis in a systematic fashion. Non-treated and treated MEFs showed distinct cell morphology, based on lipid accumulation, and expression of aP2, Lpl, and Adipoq. This finding is reflected in transcriptome data, where correlation analysis, according to distance heatmap and PCA plots, showed tight clustering between samples and corresponding replicates. The largest differences were found among treatment conditions and, to a lesser extent, with time. Thus, global changes do occur in the presence of inducers, as previously described.
Comparative pairwise analysis between treatments versus their undertreated counterpart (i.e. MEFD3plu/ MEFD3min and MEFD5plu/MEFD5min) was used to identify DEGs. A p-adj cut-off value <0.05 was applied to eliminate nonsignificant fold expression and select genes in an unbiased fashion. The approach selects significant changes regardless of fold induction. A total of 9186 and 8226 DEGs were identified in MEFD3plu/MEFD3min and MEFD5plu/MEFD5min, respectively. These slight yet variable numbers between timepoints suggest that some genes are expressed at specific times during adipogenesis. Previous studies illustrate the formation of different subtypes of adipocytes (i.e. brown, white, and beige) during adipogenesis. White adipocytes are formed earlier than brown adipocytes [2,5,62,63]. Data from the present study will require further analysis to determine if the same is true for MEFs as they differentiate into adipocytes in vitro. Moreover, linking the latter to epigenetic regulation would be useful for complete understanding of MEF adipogenic potential [11,64,65]. However, this study focused on DEGs (5314) that commonly changed globally at different timepoints. GO term analysis of these DEGs showed significant enrichment for genes involved in fat cell differentiation and adipose tissue development. The GO terms associated with fat differentiation and AT development, 119 DEGs, was further analysed. The generated heatmap showed significant upregulation exclusively in treated cells for key genes known to be involved in adipogenesis [66--66-82], including Mrap, Adig, Ffar2, Fabp4 (aP2), Slc2a4, Lrg1, Zbtb16, Rorc, Lpl, Rarres2, Wfdc21, Adrb2, Adrb3, Steap4, Retn, Adipoq, Cebpa, Cebpb, Cebpd, Lamb3, Rgs2, Scd1 and Fam57b. Interesting, when filtering key terms associated towards differentiation, there was significant enrichment terms that were linked towards fat development. This can suggest the potent effectiveness of inducing and driving the expression of pro-adipogenic genes at the biological level. However, analysis of other categories, such cellular components and/or molecular function, would further shed light mechanistically of adipogenesis in MEFs.
The magnitude of expression of DEGs might affect their capacity to affect adipogenesis. However, other factors may also play central roles in various processes, such as maintaining growth, proliferation, and functionality of adipocytes during differentiation. Hence, a PPI network of DEGs based on RRA [26] identified through the GO terms was constructed, and the top 30 genes with the highest node degrees were selected. Several ranked genes are known to have major regulatory roles in adipogenesis. These genes include protein kinase B (PKB) gene (Akt1) (the highest-ranked gene), known as a major contributor for maintaining adipocyte and adipose tissue mass through insulin signalling [83]. Moreover, glycogen synthase kinase 3-β gene (Gsk3b) was identified. This gene negatively regulates the Wnt pathway during lineage commitment towards adipogenesis in mouse embryonic stem cells. The gene product inhibits retinoic acid receptor β and negatively inhibits brown adipocyte programming [84,85]. Interestingly, cyclin D1 (Ccnd1), a cell cycle regulator that downregulates during the early phases of adipogenesis in mesenchymal stem cells [86], is downregulated specifically in treated cells in the current study. Some DEGs, significantly induced with treatment, also showed high node degrees, including Adipoq, Cebpa, Slc2a4, Fabp4 (aP2), Retn, Scd1, Lpl, Cebpb and Adrb3. These genes are known as regulators of adipogenesis and are involved in lipid metabolism, transcription, and endocrine functions in adipocytes. Consistently, genes reported in this study reflect the magnitude of signalling necessary for adipogenesis and AT as functioning tissue. Results may support the use of MEFs as a cellular model for adipocyte differentiation in vitro. However, as mentioned previously, more in-depth study of GO terms in molecular function and cellular component categories, as well as KEGG pathways, might improve understanding of adipogenesis in MEFs.
Expression changes in key genes in the DEGs list of the GO term fat cell differentiation were further analysed for expression after treatment. Cebpa and Cebpb were selected based on their identification early in landscaping adipogenesis [3,9,13,87,88] and showed expression consistent with RNA-Seq data. These transcription factors play an essential role in promoting adipogenesis, along with the master regulator Pparg [9]. However, these genes can also be detrimental to the expression of subtype specificity [88]. Such effects may be due to differential expression in a timedependent manner, as observed in this study. Specifically, some genes show variable interaction profiles with one CEBP factor depending on timepoint (day 3 vs day 5). For instance, 63 and 34 gene interactions for Cebpa were observed to be differentially expressed with MEFD3plu/ MEFD3min and MEFD5plu/MEFD5min, respectively. Nine genes were observed to be common, and 54 and 25 genes were distributed uniquely with MEFD3plu/ MEFD3min and MEFD5plu/MEFD5min conditions, respectively (Supplementary Figure 1a; Table 9). For Cebpb, 11 and 8 gene interactions were observed to be differentially expressed with MEFD3plu/MEFD3min and MEFD5plu/MEFD5min, respectively. Five genes were common, and six and three genes were distributed uniquely in MEFD3plu/MEFD3min and MEFD5plu/MEFD5min, respectively (Supplementary Figure 1b; Table 10).
The majority of genes and their interactions play pivotal roles in maintaining adipose tissue function and properties including lipid, glucose and caloric metabolism, regulation to adipocyte size and lineage commitment, epigenetic and gene transcription regulatory activities, and cell cycle regulation [86,[89][90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105] (Supplementary Figure 4). For instance, the Gsk3b gene maintains adipogenic lineage commitment during differentiation of mouse embryonic fibroblasts [84] and shows the second-highest number of degrees in the PPI network. This gene was commonly observed in the current data at both timepoints, suggesting an important and conserved role in adipocyte differentiation. The marked sequential changes in the expression of several genes during differentiation treatment suggest initiation or inhibition of various pathways for many biological functions. In light of this, several studies have noted, and demonstrated in vitro, the reduced proliferation activity during adipogenesis [86,-106-108]. For instance, the deletion of cyclin D1 (Cdk1) gene result in decrease expression of histone deacetylase (Hdac1) and increase of Pparg and therefore adipogenesis [109]. Although mitotic clonal expansion prior to adipogenesis remains to have conflicting views [110], these mechanisms have been characterized in murine models such as 3T3-L1 and mesenchymal stem cells [86,106]. In this study, similar trend of proliferation activity was observed in MEFs with Hdac1 downregulated at day 3 of differentiation along with Cdk1. These downregulation effects were restored back at day 5, suggesting that the differentiation cocktail might play a potent inhibitory role. Interestingly, the PPI data in this study showed Hdac1 to be highly ranked with 9 degrees at day 3, which also included direct interaction with Cebpa and Cebpb. Though the association of Hdac1 and Cebp family of transcription factors have been observed in MEFs [111], 3T3-L1 [112], and other cell lines [113], this study systematically identified Hdac1 having to play a pivotal role during adipogenesis in MEFs. Nonetheless, other factors identified here can attribute more to the nature aspect of adipogenesis in MEFs, which will require further analysis.

Conclusion
This study highlights roles of Cebpa and Cebpb in regulating adipogenesis through interactions of large numbers of genes. Their expression was regulated in this study in a time-dependent manner. Nonetheless, characterizing relationship among interaction of identified genes with CEBP factors (or other pioneering factors) will be essential in future work. This study has demonstrated the strength in ranking method into stressing the most influential factors driving adipogenesis in MEFs via transcriptome data. Systematic global gene expression changes during adipogenesis allow us to better understand this process for future therapeutic targeting. In addition to studying known model for adipogenesis, such as 3T3-L1 and mesenchymal stem cells, this study illustrates the use of MEFs to better understand adipogenesis from transcriptome data by identifying new key proteins that can potentially regulate adipogenesis. As such, this study has systematically noted the main factor that drive differentiation was due to reduction in proliferation activity given the downregulation of cell cycle regulated genes and its associates. Indeed, further analysis from this study data will be needed to decipher MEFs adipogenic potential.

Author contributions
Experimental designs, data analysis's, and figure preparations of study was performed by MAS. Graphic illustration of figures was performed by KA. Experiments were carried out by ME and TK. Writeup of manuscript was collectively completed by MAS, HA, MJ, and MAF.  Figure 1) highlighting common and unique genes under interactions with Cebpa (Table 9) and Cebpb (