Transcriptomic profiling of the telomerase transformed Mesenchymal stromal cells derived adipocytes in response to rosiglitazone

Differentiation of Immortalized Human Bone Marrow Mesenchymal Stromal Cells - hTERT (iMSC3) into adipocytes is in vitro model of obesity. In our earlier study, rosiglitazone enhanced adipogenesis particularly the brown adipogenesis of iMSC3. In this study, the transcriptomic profiles of iMSC3 derived adipocytes with and without rosiglitazone were analyzed through mRNA sequencing. A total of 1508 genes were differentially expressed between iMSC3 and the derived adipocytes without rosiglitazone treatment. GO and KEGG enrichment analyses revealed that rosiglitazone regulates PPAR and PI3K-Akt pathways. The constant rosiglitazone treatment enhanced the expression of Fatty Acid Binding Protein 4 (FABP4) which enriched GO terms such as fatty acid binding, lipid droplet, as well as white and brown fat cell differentiation. Moreover, the constant treatment upregulated several lipid droplets (LDs) associated proteins such as PLIN1. Rosiglitazone also activated the receptor complex PTK2B that has essential roles in beige adipocytes thermogenic program. Several uniquely expressed novel regulators of brown adipogenesis were also expressed in adipocytes derived with rosiglitazone: PRDM16, ZBTB16, HOXA4, and KLF15 in addition to other uniquely expressed genes. Rosiglitazone regulated several differentially regulated genes and non-coding RNAs that warrant further investigation about their roles in adipogenesis particularly brown adipogenesis.

approaches have many limitations which are pointing to the need for finding a new, novel, and innovative approach to treat obesity effectively [4,6]. Stem cells of different types have shown their broad capacity and effectiveness in the treatment of different diseases through their differentiation potentials. Utilizing adipose-derived stromal cells through cell-based therapy seems a promising strategy to manage obesity and related syndromes [4]. However, further understanding of adipogenesis is required for the development of effective treatment [7].
Mesenchymal Stem Cells (MSCs) are multipotent cells that has the capacity of differentiating into a variety of mesodermal cells including adipocytes [7]. Therefore, MSCs play a vital role in obesity through the generation of adipocytes, and the differentiation is considered an in vitro model of obesity [7,8]. Adipogenesis is characterized by sequential changes in the cell's gene expression profile, primarily at the transcriptional level and then differential regulation of proteins [9]. Various early, intermediate, and late markers such as mRNAs and proteins are expressed as a result of activation by several groups of transcription factors, hormones, growth factors, and extracellular matrix (ECM) proteins [9,10]. All of these modulators work in an ordered multistep process by transferring extracellular growth and differentiation signals and regulating the whole differentiation process intracellularly. MSCs will initiate to accommodate the spherical shape, enlarge and accumulate triglyceride droplets in their cytoplasm displacing the nucleus to the cell periphery, and acquire the biochemical characteristics of a mature adipocyte [11,12]. The multistep process of adipogenesis is divided into two major phases [7]. The first phase is known as the determination or the commitment phase where the multipotent MSCs commit to the adipocyte lineage and appear as pre-adipocytes. The second phase is known as terminal differentiation. Here, the pre-adipocytes are converted to mature adipocytes acquiring the full characteristics and the necessary adipocyte-specific machinery [9,13].
Adipose tissue is classically divided into two subtypes: Brown Adipose Tissue (BAT) and White Adipose Tissue (WAT) [9]. White adipocytes are the primary site of fat storage in the form of triacylglycerol in periods of energy excess, and the main fat metabolism orchestrator that works to release energy during energy deprivation [9,10]. When the energy requirements exceed the energy reserves, the stored triacylglycerol is mobilized as free fatty acids and glycerol through lipolysis [14]. Brown adipocytes, on the other hand, serve to dissipate energy through thermogenesis rather than fat storage and are relatively scarce unlike the widely distributed white adipocytes [9,10]. Given these facts, it is concluded that excess WAT is the main cause of obesity.
For effective prevention, management, and better therapeutic intervention of obesity, it is essential to study adipogenesis from progenitor cells to mature adipocytes and unravel the molecular mechanisms in such differentiation. This can be achieved by identifying the main signaling pathways and different genes that play a key role in the differentiation process. In adipose tissue, the nuclear peroxisome proliferator-activated receptor γ (PPAR-γ) is a ligand-activated transcription factor being the master regulator of BAT and WAT adipogenesis. It has vital roles in glucose and fatty acid metabolism [15]. Rosiglitazone is one of the thiazolidinediones drugs (TZDs) that was used as an antidiabetic drug and is a PPAR-γ analog [15,16]. As reported in our previous study, rosiglitazone enhanced adipogenesis by overexpression of the two transcription factors: PPAR-γ and CCAAT/enhancer binding protein α (C/EBP-α). More specifically, brown adipogenesis was enhanced by the upregulation of Early B Cell Factor 2 (EBF2) and Uncoupling protein 1 (UCP1) [17]. We reported that rosiglitazone enhances brown adipogenesis in association with the upregulation of the MAP kinase and PI3 kinase pathways. However, a deeper understanding of genes regulation during adipogenic differentiation, particularly brown adipocytes, and the effects of rosiglitazone on the transcriptomes during the differentiation is needed to be unraveled. Therefore, in this study, we investigated the transcriptomic profiles of iMSC3 and the differentiated adipocytes from iMSC3 in the presence and absence of rosiglitazone. This transcriptomic study confirmed our previous findings and further our understanding about the molecular processes that govern the adipogenic differentiation program of iMSC3, and the effects of rosiglitazone on the enhanced adipogenic differentiation, particularly brown adipogenesis.

Rosiglitazone enhances the differentiation of iMSC3 cells into adipocytes
To unravel the role of rosiglitazone in adipogenesis, the iMSC3 were differentiated in vitro into adipocytes without and with the addition of 2 μM of rosiglitazone. The morphological changes at the beginning and the end of the differentiation cycles are demonstrated in (Fig. 1). The undifferentiated iMSC3 adherent cells have fibroblast like morphology (Fig. 1A). At the end of the differentiation cycle, the adipocytes from control (Fig. 1B) and treated cells (Fig. 1C,D) were stained with Oil-O red and nile red to specifically visualize the cytoplasmic LDs formation under different experimental conditions, and DAPI to stain the nucleus. The observed morphological changes in control and rosiglitazone treated cells are characteristics of mature adipocytes. The intensity of the stain increased in adipocytes with 2 μM rosiglitazone treatment (Fig. 1C,D) in comparison with the control adipocytes (Fig. 1B) as indicated by the arrows. The stain was most intense in adipocytes derived in the presence of rosiglitazone in both induction and maintenance media (Fig. 1D). The number of lipid vesicles greatly increased and enhanced with rosiglitazone treatment. This shows that rosiglitazone enhanced adipogenesis at the morphological level. Our previous study confirmed that rosiglitazone significantly increased the lipid content of the differentiated adipocytes through lipid quantification and increased the expression of Fatty Acid Synthase (FASN) gene responsible for triglycerides synthesis [17].

mRNA sequencing, mapping and quantification
To understand the molecular mechanism of rosiglitazone in enhancing adipogenesis at the transcriptomic level, RNAseq was carried out. The sequenced mRNAs were obtained following the experimental plan depicted in (Fig. 2). To ensure the quality of downstream analysis, the sequencing raw reads were filtered to obtain clean reads by removing adaptor sequences or low-quality reads. The sequencing had effectively generated large numbers of high quality pairedend reads in all samples. All data quality is summarized in (Table S1). Spliced Transcripts Alignment to a Reference (STAR) software was used to map clean reads directly to the reference transcriptome for the differential expression gene (quantification) analysis. The summary of reads mapping to the reference genome is reported in (Table S2).

Differential gene expression analysis
The abundance of transcripts reflects gene expression level, which is calculated by the number of mapped reads and represented as Fragment Per Kilobase per Million mapped reads (FPKM) value. Read counts are proportional to gene expression level, gene length, and sequencing depth. The read counts obtained from gene expression analysis as FPKM values were used for the analysis of Differentially Expressed Genes (DEGs). The analysis was performed for each two comparison groups separately with biological replicates using the DESeq2 R package. A total of 1508 genes were found to be differentially expressed between undifferentiated iMSC3 and the fully differentiated adipocytes A vs B, among which 757 were downregulated and 751 were found to be upregulated. The genes are involved in the adipogenic differentiation of iMSC3 and consequently there is large transcriptomic changes between A and B. The comparison between the adipocytes derived in the presence of rosiglitazone added in induction media only with adipocytes derived without any addition of rosiglitazone C vs B revealed that 65 genes were downregulated and 21 upregulated giving a total of 86 DEGs. Furthermore, by comparing the transcriptomes of adipocytes derived with rosiglitazone in induction only and adipocytes derived in the presence of rosiglitazone in both induction and maintenance media C vs D, a total of 214 genes were found to be differentially expressed. Downregulated genes were 64, while the upregulated genes were 150. Surprisingly, only one significant differential expression was observed between fully rosiglitazone treated adipocytes and untreated adipocytes D vs B in FABP4 gene (Fig. 3A). Volcano plots were used to infer the overall distribution of DEGs (Fig. 3B). The top 20 DEGs in each sequenced group are listed in (Table 1). The list of all differentially regulated genes is included in (Supplementary File 2).

Co-expression analysis
The co-expression Venn diagram presents the number of genes that are both uniquely and commonly expressed within each group comparison. Comparing undifferentiated iMSC3 A and the derived adipocytes in the absence of rosiglitazone B, a total of 584 genes were found to be uniquely expressed in A and 690 genes in B sharing 12,604 genes, including many involved in the adipogenesis. When the control group B is compared to adipocytes derived in the presence of rosiglitazone added in induction media only C, the number of co-expressed genes obtained is 12,833. Notably, group B has more unique genes than C, having a total of 461 and 251 genes, respectively. On the other hand, when C transcriptome is compared to adipocytes derived in the presence of rosiglitazone added in both induction and maintenance media D, the analysis demonstrates that group D has 551 unique genes compared to 326 genes for C. Finally, the comparison of group B with D yields a total of 12,948 coexpressed genes. Group D has 361 uniquely expressed genes while B showed only 346 genes. Overall, the later pair compared showed a higher number of uniquely expressed genes (Fig. 4), in contrast to the number of DEGs within the group. The list of all uniquely expressed genes is included in (Supplementary File 3).

GO and KEGG enrichment analysis of DEGs enriched significant signaling pathways under rosiglitazone treatment
To functionally classify the differentially regulated genes and to identify their involvement in metabolic pathways, GO & KEGG enrichment analyses were performed. Through enrichment analysis of the DEGs, significant biological GO terms or pathways were found to be enriched amongst the different groups. GO and KEGG enrichment analyses were performed using ClusterProfiler software with P value < 0.05. By comparing the DEGs in undifferentiated iMSC3 with fully differentiated adipocytes A vs B, 574 significant GO terms and 6 KEGG pathways were obtained [18]. The main 20 GO terms that are enriched during the iMSC3 differentiation into adipocytes B are given in (Fig. 5A&B). According to the results, DEGs between these two groups were mainly enriched in chromosomal assembly, cell cycle pathways, and ECM related genes. The main downregulated GO terms were associated with mitosis, including mitotic sister chromatid segregation and positive regulation of cell cycle. Pointedly, KEGG enrichment analysis revealed the presence of many significant signaling pathways including PI3K-Akt (Fig. 5C).
The main GO and related terms enriched in the group C vs B were ECM related terms, cellular response to prostaglandin stimulus, and phospholipase G-protein activating protein (Fig. 6A). KEGG enrichment analysis identified regulation of PI3K-Akt pathway (Fig. 6B). This pathway is regulated with rosiglitazone treatment as reported in our previous study [17].
The main and related GO terms in C vs D enriched in the upregulated genes were lipid droplet, plasma membrane receptor complex, lipoprotein particle, protein− lipid complex, and many ECM related genes (Fig. 7A). These enriched terms do indicate that the rosiglitazone induces brown adipogenesis. KEGG pathway enrichment indicates that PPAR signaling pathway is upregulated in D due to rosiglitazone treatment (Fig. 7B).
The GO terms enriched in adipocytes derived with the addition of rosiglitazone in both induction and maintenance media D in comparison to adipocytes B were due to the expression of the FABP4. The main GO terms enriched include fatty acid binding, lipid droplet, as well

RNA-Seq validation by qRT-PCR
To confirm the differential gene expression data obtained by RNA sequencing, the expression of 13 genes were analyzed by qPCR. The genes were selected after performing DEGs analysis based on their relevance to adipogenesis. Overall, both RNA-seq and RT-qPCR showed same pattern of differential expression. The differential expression fold changes estimated by RT-qPCR for all 13 genes tested and the log transformed RNA-seq expression values (log2 fold change) were corresponding (Fig. 9).

Discussion
In this study, we investigated the changes in the transcriptome profiles of terminally differentiated adipocytes derived, with and without rosiglitazone, from the undifferentiated iMSC3 by mRNA-Sequencing. Vast transcriptomic changes were associated with the differentiation of iMSC3 into adipocytes. The rosiglitazone treatment also regulated several genes that enhanced the adipogenic differentiation, specifically brown adipocytes.
The comparison of undifferentiated iMSC3 A to adipocytes B revealed the upregulation of many adipocytes markers confirming the successful adipogenic differentiation. The list included the novel adipocytokine Retinol Binding Protein 4 (RBP4) that is known to be associated with obesity, insulin resistance, and cardiovascular diseases [19]. Several studies reported increased RBP4 expression with increased adipose tissue mass and its elevated serum level in obese human subjects [20]. Growth Differentiation Factor 15 (GDF15) is another adipokine that regulates lipid and glucose metabolism, increases insulin sensitivity, and appears to be important in maintaining body weight and preventing chronic inflammation [21,22]. Several studies showed that GDF15 levels are increased in patients with obesity and diabetes [21,23].
GO terms enrichment analysis of adipocytes B in comparison to iMSC3 revealed the upregulation of ECM related genes (Fig. 5A). The ECM is a complex network composed of different proteins, proteoglycans, and polysaccharides [24,25]. The changes in the ECM  during adipogenesis are fundamental for the morphological transformation of the cells from a fibroblastic like structure to a spherical shape [24,26]. Moreover, ECM remodeling and reorganization is essential for the regulation of adipogenesis by altering the expression of adipogenic genes, as well as for the enlargement of the existing adipocytes and the formation of new ones [26]. The ECM also provides strong external support for the mature adipocytes under strong mechanical stress; due to stored fats as triglycerides [25]. Consequently, the more mature adipocytes store fats, the more the ECM is expanded to accommodate this increase in cell volume [27]. Hence, a significant upregulation of ECM related GO terms in adipocytes B compared to undifferentiated iMSC3 cells A was observed. ECM synthesis, composition, and remodeling are structured based on the requirements for differentiation and maintaining the balance between flexibility and integrity of the tissue [28]. The downregulated GO terms in adipocytes B were chromosomal assembly and cell cycle related genes (Fig. 5B). It has been previously described that MSCs undergo a mitotic clonal expansion during early adipogenesis followed by growth arrest and adipogenic commitment [29,30]. Proliferation related genes such as cyclin D1 (CCND1), the cell cycle master regulator cyclin dependent kinase 1 (Cdk1), cell division cycle 6 (CDC6), cyclin A2 (CCNA2), polo like kinase 2 (PLK2) were clearly downregulated in differentiated adipocytes, confirming that reduced proliferative activity promotes adipogenesis as described in a previous study [30]. KEGG analysis unveiled few significant signaling pathways (Fig. 5C). PI3K-Akt signaling pathway is known to play a fundamental role in cellular processes including lipid metabolism and glucose homeostasis, protein synthesis, and cell proliferation and survival. This pathway is thought to promote lipid biosynthesis and inhibits lipolysis [31]. Evidently, this pathway was upregulated in adipocytes B. Furthermore, the transcriptome analysis of uniquely expressed genes in adipocytes B revealed many other TFs including some that are linked to adipogenesis: TFAP2E [32], POU5F1 (Oct4) [33][34][35], ZBTB16 [36,37], EGR2 [38], and MAFB [39]. These TFs along with many other genes and non-coding RNAs uniquely expressed in adipocytes warrant further investigation to find their rules in adipogenesis or obesity (Supplementary File 3). Overall, the upregulation of these markers confirms iMSC3 differentiation into adipocytes. These vast transcriptome changes in the derived adipocytes B indicate that the differentiation is associated with enormous transcriptomics changes.
To further confirm the effect of rosiglitazone treatment on the enhancement of adipogenesis, the transcriptome profiles among B, C and D were analyzed. Among the significantly enriched gene between C vs B is phosphoenolpyruvate carboxykinase 2 (PCK2) which was upregulated in C. PCK2 is a transcriptional inducer that directs the activation of phosphoenolpyruvate carboxykinase (PEPCK-C) enzyme during adipogenesis [40,41]. This enzyme catalyzes the glyceroneogenesis pathway in adipocytes that is important for fatty acid re-esterification and reduced fatty acid release and is robustly expressed in brown adipose [41,42]. As PCK2 is a peroxisome proliferator activated receptor γ response element (PPRE), it is also activated by thiazolidinediones [40]. This might explain the upregulation of PCK2 in adipocytes treated with rosiglitazone C. Within the uniquely expressed genes in adipocytes C vs B, the transcription factor PBX/knotted 1 homeobox 2 (PKNOX2) was present. This TF was reported as a novel potential regulator of browning in the bulk RNA-seq of five different mouse strains [43]. The same gene was upregulated in C compared to D. This suggests that it might have a regulatory role in the browning program that is still unexplored. In addition, the calponin 1 (CNN1) gene which was reported in beige adipocytes was also found uniquely expressed in adipocytes C [44]. This gene was later reported to be abundantly present in isolated human brown preadipocytes [45].
Adipocytes D contained FABP4, also known as adipocyte protein 2 (aP2), which is an adipogenic functional gene that is considered as an early marker of adipogenesis and is a downstream target of PPAR-γ [46]. This protein is an intracellular lipid chaperone that can reversibly bind to lipids to regulate lipid trafficking and transport to different organelles in the cell [47]. In addition, as an adipokine, FABP4 regulates glucose and lipid metabolism when released into the bloodstream, acting as a humoral factor [48]. FABP4 is known to be expressed in BATs [49]. It was shown that FABP4 can increase thermogenesis in response to both a high-fat diet and cold exposure by promoting the intracellular conversion of thyroid hormones T4 to T3 in mice BATs. Also, elevated FABP4 expression is observed in BAT of hibernating animals and cold-induced rodents [50]. The upregulation of this gene indicates more brown adipogenesis in D, Fig. 4 The Coexpression Venn diagram presenting the number of uniquely expressed genes in each group and the number of coexpressed genes in comparison with amongst the groups: undifferentiated iMSC3 A, adipocytes derived in the absence of rosiglitazone B, rosiglitazone added in induction media only C, and adipocytes derived in the presence of rosiglitazone added in both induction and maintenance media D though the gene was also expressed in C but with lower expression. FABP4 is the target of PPAR-γ and for this reason, it is affected by thiazolidinediones [47]. FABP4 delivers PPAR-γ agonists to the nucleus, thus affecting the transcription of genes that are involved in enhanced adipogenesis. It was observed previously that rosiglitazone induces increased transcription at the FABP4 locus in mouse 3 T3-L1 adipocytes [51]. As mentioned above, FABP4 is the only gene that is upregulated in D compared with B and C. This suggests that the expression of this gene is enhanced by the constant rosiglitazone treatment. This gene is so effective that it enriched several statistically significant GO terms; some of these GO terms are adipocytes and brown adipocytes related (Fig. 8). This validates our earlier finding that rosiglitazone enhances the brown adipogenesis [17].
GO and KEGG enrichment analyses confirmed rosiglitazone effect on the activation of PPAR, and regulation of PI3K-Akt signaling pathways in both treated groups C and D (Fig. 6B and Fig. 7B). The main enriched GO terms between C vs D and associate with LDs were upregulated in D (Supplementary File 4). LDs are lipid storage monolayer present in the adipocytes cytoplasm surrounded by scaffolding proteins that control lipid passage into and out of the droplets [52,53]. Perilipin 1 (PLIN1) is an LD-associated protein and is highly abundant in brown adipose tissues [54,55]. It promotes exercise-induced browning of muscle lipid, and its deficiency is observed in obese individuals [56,57]. Noteworthy, the nuclear transcription factor PPAR-γ increases the activity of PLIN1 by binding to its promoter [58]. Another upregulated gene in D is protein tyrosine kinase 2 beta (PTK2B) that regulates multiple signaling events as a member of the focal adhesion kinase (FAK) family. The increase in PTK2B protein expression was observed in cultured murine beige adipocyte differentiation. It appears that PTK2B has an essential rule in the thermogenic gene program in beige adipocytes through interacting with key adipogenic transcription factors such as PPAR-γ and C/ EBP-α [59]. Evidently, PTK2B is involved in the activation of the MAP kinase signaling pathway by stimulating JNK and ERK1/2 activity [60,61]. brown adipogenesis of iMSC3 [17]. Additionally, many ECM related GO terms were upregulated and enriched in adipocytes D (Fig. 7A). Cell-ECM interactions are reported to influence the formation of brown adipocytes and regulate their thermogenic capacity through overexpression of UCP1 [62]. Furthermore, few other studies highlighted that BAT function is regulated by its ECM [24]. This might explain the enrichment of ECM related GO terms in adipocytes D compared to adipocytes C. Therefore, a better understanding of these interactions and their role in enhancing and regulating the thermogenic capacity is needed to further explore possible therapeutic targets. This data revealed that ECM is not only involved in the differentiation of iMSC3 into adipocytes but is also involved in the enhanced differentiation of iMSC3 into adipocytes due to rosiglitazone. In addition to the DEGs, we identified Zinc Finger and BTB Domain Containing 16 (ZBTB16), Homeobox A4 (HOXA4), and Krüppel-like Factor 15 (KLF15) to be uniquely expressed in D (Supplementary File 3). ZBTB16 is known as a novel regulator of adaptive thermogenesis and is associated with the increased mitochondrial biogenesis and expression of many brown adipogenic markers [63,64]. HOXA4 was reported as a potentially important positive regulator of brown adipogenesis during the adipogenic differentiation of immortalized murine pre-adipocyte cell line [65]. Whereas KLF15 is a master regulator of adipocytes differentiation and fasting responses playing key roles in the regulation of glucose, lipids, and amino acids metabolism [66,67]. However, its physiological role in BAT needs to be unraveled [66]. Hence, the unique expression of ZBTB16 and HOXA4 do indicate that rosiglitazone enhances brown adipogenesis. Moreover, many novel noncoding transcripts were detected in rosiglitazone treated adipocytes such as lincRNAs, miRNA, snRNA, or snoRNA. A recent transcriptome study on human WATs and BATs lncRNAs revealed the roles of these non-coding RNAs in brown adipogenesis [68]. Several lncRNAs found in this study appeared to be present in our data as well, specifically in adipocytes D. The specific roles of these transcripts during adipogenesis, particularly in human adipocytes, deserve further investigations. The non-coding RNAs might affect the epigenetic changes during the differentiation that might enhance the brown adipogenesis.
Remarkably, rosiglitazone addition in both induction and maintenance media triggered the unique expression of PR domain containing 16 (PRDM16) in adipocytes D. This master transcriptional co-regulator has a crucial role in promoting the expression of brown adipocytes genes and repressing white selective genes, hence considered as one of brown adipocyte selective genes [69,70]. PRDM16 is also essential for the determination and function of beige adipocytes [69]. It was also presented that PRDM16 modulates the switch between skeletal myoblasts and brown fat cells. This regulator binds to PPAR-γ and thus stimulates brown adipogenesis [70]. It also activates the expression of thermogenic genes by co-activating PPAR-γ and PPAR-α in adipocytes [71]. Interestingly, the Trafficking Regulator of GLUT4 1 (TRARG1) was also exclusively expressed in D. This is a positive regulator of insulinstimulated GLUT4 trafficking and insulin sensitivity that works in a PI3K/Akt-dependent manner [72], suggesting the regulation of the pathways in the presence of rosiglitazone.
The upregulation of FABP4 and the unique expression of several coding and non-coding RNAs, induced by constant rosiglitazone treatment, enhanced the brown adipogenesis. These differentially and uniquely expressed genes and signaling pathways that are regulated in the presence of rosiglitazone treatment need further investigations to expand our knowledge about adipogenesis, specifically brown adipogenesis. The uniquely expressed genes and non-coding RNAs in D compared with B also suggest that rosiglitazone might regulate most of the genes at the post-translational modifications during brown adipogenesis.

Conclusion
This study provided comparisons of the transcriptomic profile of iMSC3 and adipocytes differentiated with and without rosiglitazone treatment. Moreover, it offers a reliable collection of differentially and uniquely expressed genes associated with adipogenic processes. This transcriptomic study confirmed our previous findings about the roles of rosiglitazone in the regulation of PPAR and PI3K-Akt signaling pathways during brown adipogenesis. This huge collection of data promotes broader investigations of previously studied adipogenesis and obesity related genes. Further study should be focused on the proteomics analysis to find the differentially expressed proteins that can validate our transcriptomic findings and provide further insights into the roles of rosiglitazone in the brown adipogenesis. Moreover, further study using an animal model is needed to confirm the effects of rosiglitazone on in vivo brown adipogenesis.

Nile red and DAPI staining of adipocytes
Nile red staining was performed following the protocol described by Greenspan et al. with modifications [78]. In brief, nile red stock solution of 1 mg/mL concentration was obtained by dissolving 5 mg of nile red powder (Sigma-Aldrich N3013, USA) in 5 mL of acetone. The stock solution was diluted to 1:100 nile red staining solution in 1 mM trizma-maleate (Sigma-Aldrich T3128, USA) and 3% w/v Polyvinylpyrrolidone (Sigma-Aldrich P2307, USA). The differentiated adipocytes were washed with PBS (Sigma-Aldrich D8537, USA) and fixed with 4% paraformaldehyde for 1 h. Then, the cells were stained with nile red to visualize the lipid droplets and DAPI (Life Technologies P36930, USA) to stain the nucleus. Olympus fluorescent microscope was used to observe the stained cells and imaged using cellSens Standard software.

Oil-O red staining of adipocytes
Oil-O red staining was performed following Aguena et al. protocol with modifications [79]. The differentiated adipocytes were washed with PBS, fixed with 4% paraformaldehyde for 1 h, and incubated for 15 min at room temperature with 60% isopropanol solution. Then, the cells were air-dried, stained with Oil-O red staining solution (Sigma-Aldrich O1391, USA) (60% Oil-O Red stock solution and 40% distilled water) for 1 h, and washed with distilled water to remove excess stain. The stained cells were visualized under inverted microscope using Optika Vision Lite software.

RNA extraction and quantification
Total RNAs were extracted from all experimental groups A, B, C, and D. The harvested cells were first lysed using Qiazol lysis reagent. Then, the extraction was carried out using miRNeasy extraction kit following the manufacturer's instructions (Qiagen 217,004, Germany). The RNA quantity and purity were determined using Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, USA).

RNA quality and integrity assessment
To guarantee the reliability of the data, quality control (QC) was performed at each step. RNA degradation and contamination were assessed on 1% agarose gel. Then, RNA purity was checked using NanoPhotometer spectrophotometer (Implen, USA). To obtain a quantitative assessment of RNA integrity and evaluate the RNA quantity, Agilent BioAnalyzer 2100 system was used. The samples concentration was first unified to 500 ng/ul, then proceeded with the Agilent RNA 6000 Nano Kit protocol as per the manufacturer's instructions (Agilent Technologies 5067-1511, USA). The degradation was scored for each sample and represented as RNA Integrity Number (RIN) value. Only samples with RIN > 5.0 were used for subsequent library construction (Table S3).

Library preparation for Transcriptome sequencing
A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations. NEB-Next ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, USA) was used to generate sequencing libraries following manufacturer's recommendations. To attribute sequences to each sample, index codes were added. Briefly, poly-T oligo-attached magnetic beads were used to purify mRNA from total RNA. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). In order to select cDNA fragments, preferably of 150-200 bp in length, the library fragments were purified using AMPure XP system (Beckman Coulter, USA). Then, 3 μl of USER Enzyme (NEB, USA) was used with size-selected, adaptor ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C. After that, PCR was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Finally, the PCR products were purified (AMPure XP system) and library quality was assessed using the Agilent Bioanalyzer 2100 system.

Clustering and sequencing
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using PE Cluster Kit cBot-HS (Illumina, USA) according to the manufacturer's instructions. After cluster generation, the library was sequenced using Illumina NovaSeq 6000 platform (Illumina, USA) and pairedend reads were generated in 300 cycles. All generated sequencing data are deposited in Gene Expression Omnibus (GEO) database with the accession number (GSE171826).

Data analysis Quality control
To ensure the quality and reliability of data analysis, the raw reads obtained in FASTQ format were first processed through fastp. Reads containing adapter sequences and poly-N sequences were removed from the raw data along with low-quality reads.
Simultaneously, Q20, Q30 and GC content of the clean data were calculated. All the downstream analyses were performed based on high quality clean data.

Mapping to reference genome
Reference genome and gene model annotation files were directly downloaded from genome website browser (NCBI, UCSC, and Ensembl). Using STAR software, the paired-end clean reads were aligned to the reference genome. This software is based on a previously undescribed RNA-seq alignment algorithm that uses sequential maximum mappable seed search in uncompressed suffix arrays followed by seed clustering and stitching procedure. Compared to other RNA-seq aligners, STAR exhibits better alignment precision and sensitivity for both experimental and simulated data [80]

Quantification
The gene expression level was estimated by the abundance of transcript mapped to genome or exon. Fea-tureCounts was used to count the reads number mapped to each gene. The FPKM value was calculated for each gene based on the length and reads count mapped to this gene.

RT-PCR Primer Seq (Group C vs D) References
insulin like growth factor binding protein 2 (IGFBP2) GCC CTC TGG AGC ACC TCT ACT CAT CTT GCA CTG TTT GAG GTT GTA C [90] ras related dexamethasone induced 1 (RASD1) CCA CCG CAA GTT CTA CTC CAT CCA GGA TGA AAA CGT CTC CTGT [91] fatty acid binding protein 4 (FABP4) GCC AGG AAT TTG ACG AAG TCAC TTC TGC ACA TGT ACC AGG ACAC [92] prostaglandin I2 synthase (PTGIS) CTG GTT GGG GTA TGC CTT GG TCA TCA CTG GGG CTG TAA TGT [93] apolipoprotein L3 (APOL3) GCA AGG GAC ATG ATG CCA GA AAG AGT TTC CCC AAG TCA AGAGG [94] phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 beta (PIK3C2B) CAG GCT TCA AGA GGC ACT CA TGG TCA TCA TTC ACC GTC CG [95] GAPDH AGG GCT GCT TTT AAC TCT GGT CCC CAC TTG ATT TTG GAG GGA [17] having three biological replicates, was performed using DESeq2 R package. DESeq2 provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution [81]. Using the Benjamini and Hochberg's approach for controlling the False Discovery Rate (FDR), the resulting P values were adjusted. Genes with an adjusted P value < 0.05 found by DESeq2 were considered as differentially expressed for the three biological replica.

GO and KEGG enrichment analysis
Enrichment analysis enables us to attribute biological functions or pathways that are significantly associated with the DEGs. GO enrichment analysis of DEGs was implemented by the clus-terProfiler R package. GO terms with corrected P value < 0.05 were considered significantly enriched by differential expressed genes. R package clusterProfiler was used to test the statistical enrichment of differential expression genes in KEGG pathways. Those terms with adjusted P value < 0.05 were considered as significantly enriched.

cDNA synthesis and RNA-Seq validation by quantitative RT-PCR
cDNA was synthesized from the total RNA using QuantiTect Reverse Transcription kit following the manufacturer's protocol (Qiagen 205,311, Germany). To confirm the transcriptome results, a total of 13 DEGs were selected and analyzed by qPCR using the same samples used for the RNA sequencing. The qPCR primer sequences for target genes are presented in (Table 2). QuantStudio3 Real-Time PCR System (Applied Biosystems, USA) was used to perform realtime gene expression study and using Maxima SYBR Green/ROX qPCR Master Mix (2x) (Thermo Fisher Scientific K0223, USA). GAPDH was used to normalize the obtained expression levels. The results were analyzed using Design and Analysis Software v1.5.1 (Thermo Fisher Scientific, USA), and the relative expression of genes was calculated using the 2 −ΔΔCT method [82].

Statistical analysis
All experiments were performed in triplicates and results were expressed as the mean ± standard deviation. Statistical significance was analyzed using Graph-Pad Prism 9 software to perform unpaired two-tailed student t-test considering significance at P value < 0.05.